AD-A092  bl3 


E-SYSTEMS  INC  DALLAS  TEX  GARLAND  DIV 
STUDY  OF  DIGITAL  MATCHING  OF  DISSIMILAR  IMAGES. <U> 


UNCLASSIFIED 


OCT  80  D  S  RHINES 


DAAK70-79-C-0 


DBG  file  copi 


ETL-0244 


'V 


ADA092813 


FINAL  REPORT,  STUDY  OF  DIGITAL  MATCHING  OF  DISSIMILAR  IMAGES 


Don  S.  Rhines 
E-SYSTEMS,  Inc. 
Garland  Division 
P.O.  Box  226118 
Dallas,  Texas  75266 


31  October  1980 


Approved  for  public  release;  distribution  unlimited 


Prepared  for 

U.S.  Army  Engineer  Topographic  Laboratories 
Fort  Belvoir,  Virginia  22060 


c 

£ 

W  v, 
iiSQ  ’  i  ' 


A 


80  1201  163 


Destroy  this  report  when  no  longer  needed. 
Do  not  return  it  to  the  originator. 


The  findings  in  this  report  are  not  to  be  construed  as  an 
official  Department  of  the  Army  position  unless  so  designated 
by  other  authorized  documents. 


The  citation  in  this  report  of  trade  names  of  commercially 
available  products  does  not  constitute  official  endorsement  or 
approval  of  the  use  of  such  products. 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Wh,n  Data  Bnfnd) 


■  j  REPORT  DOCUMENTATION  PAGE 

READ  INSTRUCTIONS 

BEFORE  COMPLETING  FORM 

1.  REPORT^imOEN  2.  GOVT  ACCESSION  NO. 

ETL-0244  •  Ah- A  Mi, 

9.  RECIPIENT'S  CATALOG  NUMBER 

Ul 

4.  TITLE  (and  Subtlti,)  {(,' l~ - - -  ' 

Final  Report,  Study  of  Digital  Matching 
of  Dissimilar  Images  * 

s.  type  of  report  a  period  cov^«ed 

Contract  Report 

S.  PERFORMING  ORG.  REPORT  NUMBER 

7.  AUTHOR(a) 

)  Don  S./Rhines  ? )  /  /  •  j-i»*  ! 

_ 1  . .  '  y 

8.  CONTRACT  or  GRANT  NUMBER^*; 

)  DAAK70-79-C-0235  /  ~ 

9.  PERFORMING  ORGANIZATION  NAME  AND  ADDRESS 

E-SYSTEMS,  Inc.,  Garland  Division 

P.0.  Box  226118  v 

Dallas.  Texas  75266 

10.  PROGRAM  ELEMENT.  PROJECT.  TASK 
AREA  ft  WORK  UNIT  NUMBERS 

iiyji  ■ 

It.  CONTROLLING  OFFICE  NAME  AND  ADDRESS 

U .  S .  Army  Engineer  Topographic  // J 

Laboratories,  Fort  Belvoir,  Virginia  7 

22060  s 

12.  REPORT  DATE 

31  October  1980 

13.  NUMBER  OF  PAGES _ ' 

98 

14.  MONITORING  AGENCY  NAME  S  AOORESS<7(  dJM.rant  from  Controlling  Otlic,) 

IS.  SECURITY  CLASS,  (ot  tht»  report) 

UNCLASSIFIED 

is*,  declassification/ DOWNGRADING 
SCHEDULE 

16.  DISTRIBUTION  STATEMENT  (of  thle  Report) 

Approved  for  public  release;  distribution  unlimited 

17.  DISTRIBUTION  STATEMENT  (of  the  ebetrect  entered  in  Block  30,  If  different  from  Report) 

DTIC 

f|ELECT,:  r--\ 

18.  SUPPLEMENTARY  NOTES  u  isou 

A 

18.  KEY  WOROS  ( Continue  on  reverie  eide  If  neceeeery  end  Identify  by  block  number) 

Matching,  Dissimilar  Images,  Physical  Commonalities,  Feature 
Extraction,  Similarity  Measures 

20.  ABSTRACT  fCaoUtaua  ea  revere*  ft  neee emery  end  Ideality  by  block  number) 

This  final  report  presents  the  results  of  a  study  con¬ 
ducted  for  the  U.S.  Army  Engineering  Topographic  Laboratories 
on  the  digital  matching  of  dissimilar  images.  This  report 
develops  a  practical  approach  for  the  digital  determination 
of  corresponding  points  on  dissimilar  images.  This  approach 
could  be  used  to  register  a  large  number  of  points  automatically v 

DO  I  jam  T»  1473  EDITION  OF  »  MOV  CS  is  OBSOLETE 


UNCLASSIFIED 


SECURITY  CLASSIFICATION  OF  THIS  PAGE  (Nh«n  Daia  BmaraV) 


3 


CLASSIFIED  _ 


SECURITY  CLASSIFICATION  OF  THIS  PAGElTOian  Dmlt  Knfrnd) 


SECURITY  CLASSIFICATION  of  this  FAOEf»?i»n  D»t*  En<tr»d) 


1 


PREFACE 


This  report  is  generated  under  Contract  DAAK70-79-C-0235 
for  the  U.S.  Army  Engineer  Topographic  Laboratories,  Fort 
Belvoir,  Virginia  22060  by  E-SYSTEMS,  Inc.,  Garland  Division, 
Dallas,  Texas  and  submitted  as  ETL-0244.  The  Contract 
Officer's  Representative  was  Mr.  Tom  Blackburn. 


5 


TABLE  OF  CONTENTS 

PAGE 

1 . 0  INTRODUCTION  11 

2.0  GEOMETRIC  CORRECTIONS  17 

3.0  COMPUTER  ASSISTED  MANUAL  IMAGE  REGISTRATION  23 

4.0  TEXTURE  MATCHING  CONCEPT  DEMONSTRATION  35 

4.1  Conventional  Correlation  Experiments  36 

4.1.1  No  Modification  Example  40 

4.1.2  Four  Modification  Examples  42 

4.2  The  Texture  Matching  Approach  51 

4.2.1  Texture  Matching  Example  57 

4.2.2  Resolution  Factors  61 

4.3  Some  General  Conclusions  66 

5.0  TEXTURE  MATCHING  REFINEMENTS  69 

5.1  Some  Simple  Transformations  70 

5.1.1  High  Pass  Filtering  70 

5.1.2  Edge  Enhancement  Operation  71 

5.2  Texture  Measures  from  Spectral  Analysis  71 

5.2.1  Four-by-Four  Pixel  Neighborhood  76 

5.2.2  Eight-by-Eight  Pixel  Neighborhood  77 

5.2.3  Low  Frequency  Background  Subtraction  83 

5.3  Formal  Definition,  Texture  Matching  65 

5.3.1  Single  Texture  Measure  85 

5.3.2  Multiple  Texture  Measures  87 

5.4  Miscellaneous  Corrections  88 

5.4.1  Noise  Characteristics  88 

5.4.2  Spatial  Acceptance  89 

5.4.3  Intensity  Distributions  90 

5.4.4  Coherency  Effects  92 

5.5  Match  Point  Accuracy  94 

6.0  SUMMARY  95 


7 


LIST  OF  TABLES 


Table  No.  Title 

4-1  Summary,  Match  Point  Errors,  Texture 

Matching  Examples 

5_I  Texture  Measure  Components,  8X8  Neighborhood 

Case 


Page 

65 


79 


8 


* 


LIST  OF  FIGURES 

Figure  No.  Title  Page 

4-1  Sample  Aerial  Image  37/38 

4-2  Correlation  Example,  No  Modification  Added  41 

4-3  Correlation  Example,  Noise  Modification  Added  43 

4-4  Correlation  Example,  Low  Pass  Filter  45 

Modification  Added 

4-5  Nonlinear  Intensity  Transformation  47 

4-6  Correlation  Example,  Nonlinear  Intensity  49 

Transformation  Modification  Added 

4-7  Correlation  Example,  Edge  Enhancement  52 

Distortion  Added 

4-8  Error  in  Correlation  Peak  Position,  Edge  53/54 

Enhancement  Modification  Added 

4-9  Sample  Aerial  Image  Modified  by  an  Edge  55/56 

Enhancement  Operation 

4-10  Texture  Matching  Example,  One-by-One  Pixel  59 

Replacement  Resolution 

4-11  Texture  Matching  Example,  Failure  Case  60 

4-12  Texture  Matching  Example,  Two-by-Two  Pixel  63 

Replacement  Resolution 

4- 13  Texture  Matching  Example,  Four-by-Four  Pixel  64 

Replacement  Resolution 

5- 1  Division  of  Frequency  Components  Defining  the  78 

Low,  Medium,  and  High  Frequency  Texture 
Measures 


9 


LIST  OF  FIGURES  (Continued) 


Title  Page 

Division  of  Frequency  Components  Defining  the  80 

X,  XY,  and  Y  Texture  Measures 

Definition  of  an  Independent  Set  of  Frequency  82 

Components 

Low  Frequency  Background  Subtraction  Example  84 

Example,  Matching  Spatial  Acceptance  with  91 

Low  Pass  Filtering 

Typical  Intensity  Distributions  93 


10 


1.0 


Introduction 


This  final  report  presents  the  results  of  a  study 
conducted  for  the  U.  S.  Army  Engineering  Topographic  Laboratories 
on  the  digital  matching  of  dissimilar  images.  This  report  develops 
a  practical  approach  for  the  digital  determination  of  corresponding 
points  on  dissimilar  images.  This  approach  could  be  used  to  register 
a  large  number  of  points  automatically  in  a  reasonably  short 
period  of  time.  The  algorithms  required  are  presented  in  a  manner 
that  can  be  coded  in  FORTRAN  IV  and  tested  on  the  DIAL  facility 
at  USAETL. 

For  similar  images,  conventional  correlation  techniques 
can  be  used  to  digitally  match  images.  With  the  measurement  of 
many  corresponding  points,  contour  information  can  be  developed 
using  this  approach.  However,  for  the  technique  to  work  with  any 
degree  of  reliability,  requires  that  special  optical  cameras 
be  used  and  that  the  taking  conditions  be  precisely  controlled. 

The  photographs  must  be  taken  at  the  same  time  of  day  from  just 
the  right  perspectives. 

This  report  deals  with  the  study  of  techniques  that 

would  permit  digital  matching  between  a  much  broader  class  of 

images.  It  is  desired  to  be  able  to  match  images  from  different 

detector  sources.  Optical,  infrared,  and  radar  detectors  are 

examples  of  the  bro^d  source  of  images  it  is  desired  to  be  able 

to  match.  Similarly,  it  is  desired  to  be  able  to  match  images 

from  a  single  sensor  when  the  images  are  taken  from  very  different 

perspectives  or  under  different  illumination  conditions.  The  most 

difficult  class  of  images  to  match  are  radar  images.  This  is  because 

of  the  coherent  nature  of  the  imaging  techniques.  Several  specific 

suggestions  are  included  in  this  report  to  help  deal  with  radar 
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images . 


There  are  two  different  approaches  that  can  be  taken  to 
develop  techniques  to  match  dissimilar  images.  One  could  attempt 
to  use  a  pattern  recognition  approach  or  a  texture  matching  approach. 

In  a  pattern  recognition  approach,  items  in  the  scene  are  identified 
through  feature  extraction.  In  a  texture  matching  approach,  it  is 
assumed  that  between  matching  areas  of  the  images  there  will 
remain  matching  textural  components.  With  the  texture  matching 
approach,  the  difficult  task  of  item  identifiction  is  unnecessary. 

In  this  study  the  texture  matching  approach  to  matching  dissimilar 
images  has  been  used. 

In  most  situations  there  will  exist  some  information  about 
the  taking  conditions  used  to  generate  the  images  of  interest. 

This  information  could  include  the  approximate  position  of  the 
camera  including  altitude  and  the  aiming  direction  used.  To  reduce 
the  dimensions  of  the  matching  problem,  it  is  strongly  suggested 
that  whatever  a  priori  information  is  available  be  used  initially 
in  the  matching  process  to  generate  geometric  corrections.  In 
Section  2  a  general  technique  is  presented  to  develop  geometric 
corrections.  This  technique  can  be  applied  to  optical  or  radar 
images,  and  can  include  the  use  of  contour  information  if  available. 
After  the  application  of  these  corrections,  the  dominant  misregistration 
components  will  be  translational.  Scale  and  rotational  differences 
will  normally  be  small.  Both  global  and  local  translational 
differences  will  exist.  The  global  or  overall  translational 
differences  come  from  the  residual  uncertainties  in  aim  point 
determination.  Local  translational  differences  are  due  to  the 
combined  effects  of  the  perspectives  used  and  uncorrected  terrain 
effects.  It  is  generally  assumed  in  developing  matching  algorithms 

that  geometric  corrections  of  this  type  have  been  applied. 
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There  are  two  extreme  applications  to  which  the  matching 
of  dissimilar  images  could  be  used.  On  the  one  hand,  it  could  be 
desired  to  match  thousands  of  points  for  the  development  of  a 
contour  data  base.  On  the  other  hand,  it  could  be  desired  to 
measure  a  few  match  points  from  each  image  in  order  to  approximately 
register  the  images.  After  approximate  registration,  the  global 
translational  differences  between  the  two  images  have  been  removed, 
but  the  local  translational  differences  would  still  remain.  The 
approximate  registration  of  images  might  be  required  in  a  map 
updating  application.  Also,  approximate  registration  of  the  images 
could  be  used  as  the  first  step  of  a  contour  measuring  process. 

The  texture  matching  approach  that  has  been  developed  in  this 
report  could  be  used  in  either  of  these  extreme  applications. 

In  addition  to  the  texture  matching  approach,  another 
approach  is  developed  for  use  in  the  approximate  registration  of 
images.  This  approach  is  described  in  Section  3,  and  is  a 
semi-automatic  approach  which  makes  use  of  the  human  ability  to 
find  match  points  between  differing  images.  The  approach 
provides  computer  and  special-purpose  hardware  assistance  to  the 
operator  who  is  required  to  measure  matching  points  from  each  of 
the  images.  This  approach  is  referred  to  as  computer  assisted 
manual  image  registration  (CAMIR) .  This  system  could  be  readily 
constructed  given  the  current  .state-of-the-art  of  digital 
hardware . 

The  texture  matching  approach  has  the  advantage  of 

being  completely  automatic.  The  basic  concept  is  to  replace  the 

actual  sample  values,  which  are  very  dependent  upon  the  collection 

system,  with  derived  values  which  are  independent  of  the  collector 

system.  With  system  independence  and  with  the  application  of  gross 
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geometrical  corrections,  image  pairs  can  be  registered  by  con¬ 
ventional  correlation  techniques.  The  key  to  this  approach  is 
the  development  of  specific  conversion  algorithms  for  the  pixels 
of  each  system  of  interest,  so  that  the  resultant  pixel  values 
are  effectively  system  independent.  Although  the  general  procedure 
can  be  readily  described  and  implemented,  and  although  there  is 
considerable  experience  with  specific  transformations,  the  actual 
pixel  replacment  algorithms  will  have  to  be  empirically  tailored 
for  specific  collection  systems. 

In  Section  4  the  texture  matching  concept  is  demonstrated 
A  single  aerial  image  is  used.  It  is  found  that  when  this  image 
is  modified  with  the  use  of  an  edge  enhancement  operation,  match 
points  can  no  longer  be  obtained  using  a  conventional  correlator. 
The  pixels  from  both  the  original  and  modified  images  are 
replaced  with  a  measure  of  local  texture.  The  texture  measure 
used,  referred  to  as  the  local  busyness  measure,  uses  the  local 
AC  energy  present  in  the  neighborhood  of  the  replaced  pixel.  It 
is  then  found  that  conventional  correlation  techniques  can  be  used 
to  accurately  match  the  images.  Factors  affecting  match  point 
resolution  are  further  discussed  in  this  section.  It  is  found  that 
it  is  not  necessary  to  replace  each  and  every  pixel  with  a  texture 
measure  in  order  to  maintain  match  point  resolution.  The  frequency 
with  which  pixel  intensity  values  are  replaced  is  referred  to  as 
the  pixel  replacement  resolution. 

The  texture  matching  approach  is  further  developed  in 
Section  5.  Derivation  of  local  texture  measures  through  the  use  of 
local  spectral  analysis  is  discussed.  The  local  spectral  analysis 
is  performed  with  the  use  of  the  two  dimensional  discrete  Fourier 
transform.  This  is  a  powerful  approach  which  Can  be  used  to 


generate  a  broad  class  of  texture  measures.  The  texture  matching 
approach  is  formally  defined,  and  the  technique  is  generalized  to 
permit  the  use  of  multiple  texture  measures.  Finally,  a  variety 
of  special  problems  which  can  interfere  with  texture  matching  are 
discussed.  Methods  are  discussed  for  dealing  with  data  noise 
characteristics,  different  detector  spatial  acceptances,  drastically 
different  intensity  distributions,  and  coherency  effects  found  in 
dealing  with  imaging  radars. 
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GEOMETRIC  CORRECTIONS 

In  our  approach  to  the  problem  of  matching  a  large 
number  of  points  between  the  dissimilar  images,  the  first  step 
is  to  make  all  possible  geometric  corrections  based  upon  a  priori 
knowledge.  In  this  section  geometric  correction  fundamentals  are 
outlined. 

The  first  step  is  to  derive  a  relationship  between  posi¬ 
tions  in  the  f ield-of-view  of  the  taking  system  and  the  correspond¬ 
ing  points  onto  which  they  are  projected  by  the  taking  system.  Let 
positions  in  the  f ield-of-view  of  the  system  or  in  object  space  be 
denoted  by  the  coordinates  x,  y,  and  z.  Let  positions  on  the  pro¬ 
jection  surface  be  denoted  by  the  coordinates  a  and  B.  Then  the 
desired  relationships  are  of  the  form: 


fa(x,y,z) 

(2-1) 

f g (x,y,z) 

(2-2) 

To  derive  these  relationships,  the  taking  parameters  which  would 
need  to  be  known  would  include  the  position  and  altitude  of  the 
system,  the  aiming  direction  of  the  system,  as  well  as  other 
specifics  of  the  projection  process  used  by  the  system.  In  some 
cases  a  simple  analytic  expression  can  be  fairly  easily  derived. 
An  example  of  this  would  be  a  snapshot  camera  with  no  lens  dis¬ 
tortions  present.  In  other  cases  much  more  elaborate  procedures 
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would  be  necessary  to  derive  the  appropriate  relationship.  An 
example  of  this  case  would  be  Landsat  satellite  imagery.  In 
this  case  the  satellite  flight  path  and  motion  used  during  the 
exposure  must  be  taken  into  account. 

The  next  step  is  to  determine  a  grid  of  points  on  a 
surface  in  object  space.  The  original  scene  is  distorted  due 
to  the  taking  perspective  and  characteristics  of  the  taking  sys¬ 
tem.  We  would  like  to  determine  the  appropriate  pixel  intensi¬ 
ties  in  an  undistorted  space,  and  the  grid  of  points  required 
gives  the  positions  at  which  we  would  like  to  know  the  pixel 
intensities.  Let  position  on  this  object  surface  be  denoted  by 
the  coordinates  r  and  s.  Since  we  are  developing  information 
for  use  in  map  making,  the  coordinates  r  and  s  will  usually  cor¬ 
respond  with  appropriate  map  projection  coordinates.  An  example 
would  be  the  position  on  the  surface  of  the  earth  as  specified 
by  longitude  and  latitude.  If  detailed  contour  information  is 
available,  then  the  surface  to  be  used  can  correspond  with  the 
actual  surface  of  the  earth.  If  detailed  contour  information 
is  not  known,  then  an  approximation  to  the  earth's  surface  at 
the  appropriate  average  altitude  can  be  used.  If  the  area  covered 
is  small  (less  than  a  square  mile) ,  then  a  flat  plane  would  be 
adequate  for  this  surface.  For  larger  area  coverage,  then  a 
sphere  or  oblate  spheroid  surface  would  be  an  appropriate  ap¬ 
proximation.  For  the  choice  of  grid  points  made,  then  for  each 
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grid  point  its  position  in  object  space  must  be  determined. 
In  functional  form,  the  required  relationship  is  given  by: 


X  = 

f  (r, s) 

(2-3) 

X 

y  = 

fy(r,s) 

(2-4) 

z  = 

fz(r,s) 

(2-5) 

The  next  step 

is  to  determine  for  each 

grid  point  the 

corresponding  position 

in  the  projection  plane. 

Substituting 

Equations  (2-3) ,  (2-4) , 

and  (2-5)  into  (2-1)  and 

(2-2)  ,  we 

obtain  the  desired  relationships : 

a  = 

fa(r,s) 

(2-6) 

5  = 

£8(t,sl 

(2-7) 

For  many  cases  it  should  be  possible  to  derive  this  relationship 
directly.  This  can  be  done  as  long  as  the  object  surface  can 
be  described  by  some  simple  functional  form.  The  one  clear  ex¬ 
ception  is  the  case  in  which  contour  information  is  being  used. 

In  this  case,  the  altitude  information  would  be  contained  in  a 
grid  of  altitude  values.  In  this  case  it  would  not  be  reasonable 


to  derive  explicitly  relationships  (2-6)  and  (2-7) .  Instead,  a  two 
step  procedure  would  be  required.  First,  for  each  grid  point, 
its  position  in  object  space  must  be  determined  by  interpolation 
from  the  altitude  data  base.  Next,  the  corresponding  position 
on  the  projection  surface  can  be  determined  using  Equations  (2-1) 
and  (2-2) . 

In  the  final  step,  the  intensity  values  associated  with 
each  of  the  grid  points  which  have  been  mapped  onto  the  projection 
plane  are  determined  in  an  interpolation  process.  The  interpola¬ 
tor  requires  the  position  of  the  grid  point  with  respect  to  the 
pixel  locations  and  the  value  of  the  intensity  levels  associated 
with  the  pixel  locations  in  the  immediate  neighborhood  of  the 
projected  grid  point.  In  many  applications  a  2x2  point  linear 
interpolator  should  give  adequate  results.  Performance  can  be 
improved  by  going  to  a  higher  order  interpolator  using  specially 
developed  coefficients.  At  E-Systems,  we  have  extensive  experi¬ 
ence  in  developing  coefficients  for  higher  order  interpolators 
should  high  performance  be  required. 

When  hardware  implementations  of  geometric  corrections 
are  considered,  it  becomes  clear  that  including  contour  informa¬ 
tion  adds  considerably  to  the  complexity  of  the  problem.  If 
contour  information  is  to  be  used,  the  process  probably  should 
be  divided  up  into  a  two-stage  process.  Initially,  the  geometric 
corrections  should  be  generated  without  the  use  of  the  contour 


20 


information.  The  image  can  then  be  registered  with  respect  to 
a  reference  image.  Using  this  process  the  knowledge  of  the  taking 
parameters  can  be  refined.  Without  this  initial  step,  knowledge 
of  the  taking  parameters  would  not  normally  be  sufficiently  ac¬ 
curate  to  warrant  the  use  of  contour  information.  In  the  second 
stage,  the  full  geometric  corrections  with  the  use  of  the  contour 
information  can  be  generated. 

There  is  one  other  special  problem  that  can  arise  with 
the  use  of  contour  information.  In  general,  if  the  surface  is 
bumpy  or  mountainous  and  the  look  angle  well  off  nadir  towards 
the  horizon,  then  there  will  be  regions  on  the  object  surface 
which  will  be  hidden  from  the  detector.  For  geometrically  correct¬ 
ed  optical  imagery,  these  regions  should  be  flagged  in  some  un¬ 
ambiguous  manner. 

Consideration  of  the  addition  of  contour  information 
should  depend  on  the  class  of  images  to  be  used  in  making  maps. 

In  many  map  updating  operations,  contour  information  should  be 
frequently  available  and  be  useful.  The  more  mountainous  the 
terrain,  the  larger  will  be  the  distortions  that  can  be  generated. 
In  general,  the  more  an  optical  detector  looks  off  nadir  toward 
the  horizon,  the  more  severe  will  be  these  distortions. 
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3.0  COMPUTER  ASSISTED  MANUAL  IMAGE  REGISTRATION 

To  solve  the  problem  of  matching  dissimilar  images,  it 
is  suggested  that  the  first  step  is  to  bring  images  into  approxi¬ 
mate  registration.  In  this  section  we  describe  a  technique  that 
can  be  used  to  obtain  approximate  registration.  This  technique 

makes  use  of  high-speed  digital  approaches  to  applying  geometric 
corrections,  and  makes  use  of  the  exceptional  ability  of  a  human 
operator  to  quickly  solve  complex  pattern  recognition  problems. 
This  approach  is  referred  to  as  computer  assisted  manual  image 
registration  (CAMIR) .  The  state-of-the-art  of  digital  hardware 
has  advanced  to  such  a  point  that  this  system  could  be  readily 
constructed.  It  is  felt  that  this  approach  offers  the  most 
practical  immediate  solution  to  the  problem  of  quickly  bringing 
into  approximate  registration  dissimilar  images. 

First,  we  will  describe  the  operation  and  performance 
capabilities  of  the  system  that  is  envisioned.  The  system  will 
require  both  the  use  of  a  general-purpose  computer  and  the  use 
of  high-speed  special-purpose  hardware.  The  particular  approach 
required  to  correct  geometric  distortions  will  be  described. 

This  unique  approach  permits  the  real-time  implementation  of  the 
geometric  corrections  for  any  type  of  detector  using  high-speed 
digital  hardware.  Details  of  the  hardware  are  beyond  the  scope 
of  this  study,  but  cursory  analysis  indicates  that  the  hardware 
requirements  are  well  within  the  current  state-of-the-art. 
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To  view  each  of  the  two  dissimilar  images,  the  operator 
will  view  a  CRT  display.  The  CRT  would  normally  display  just 
a  small  portion  of  the  image.  The  operator  would  have  complete 
freedom  to  explore  and  examine  a  given  image.  This  freedom  would 
include  the  ability  to  translate,  to  rotate,  and  to  both  magnify 
or  demagnify  the  sub-image  being  viewed.  Updating  of  the  sub¬ 
image  would  be  performed  at  a  sufficiently  high  rate  to  assure 
smooth  presentation  of  manipulations  to  the  display.  Geometric 
corrections  based  on  a  priori  knowledge  are  normally  applied 
before  the  image  is  projected  onto  the  CRT. 

To  approximately  register  the  two  images,  the  operator 
would  select  three  or  more  corresponding  points.  For  best  per¬ 
formance  these  points  should  be  selected  as  much  as  possible  at 
the  extremes  of  the  images,  as  widely  separated  as  possible.  In 
order  to  make  the  measurements,  the  operator  would  bring  a  cursor 
projected  on  top  of  the  CRT  display  into  coincidence  with  the 
point  to  be  measured. 

After  three  corresponding  points  in  each  image  have 
been  measured,  the  geometric  corrections  would  be  refined  to 
bring  the  two  images  into  approximate  registration.  The  three 
points  which  have  been  measured  define  a  plane.  All  points  which 
lie  in  this  plane  will  be  now  exactly  registered.  Points  lying 
above  or  below  this  plane  would  be  shifted  with  respect  to  one 
another  from  one  image  to  the  next.  The  amount  of  shift  would 
depend  on  how  far  the  points  are  from  the  registration  plane  and 
on  the  difference  in  the  taking  perspectives.  For  many  pairs  of 

24 


-  r 


images  the  remaining  shifts  present  over  most  of  the  field-of- 
view  would  be  quite  small.  This  fact  can  now  be  used  in  develop¬ 
ing  further  techniques  to  register  many  points  between  the  images. 

If  it  were  desirable,  it  is  possible  at  this  point  to 
develop  a  system  which  would  simultaneously  project  the  two  images 
for  stereo  viewing.  As  long  as  the  taking  perspectives  are  ap¬ 
propriate,  the  two  images  could  be  from  completely  different 
taking  systems.  For  stereo  viewing  of  any  portion  of  the  image, 
no  further  manually  applied  corrections  would  be  necessary. 

There  are  several  other  highly  desirable  features  which 
could  be  included  under  the  real-time  control  of  the  operator. 

One  of  these  is  the  ability  to  adjust  the  correction  being  ap¬ 
plied  to  account  for  the  spatial  frequency  response  found  in 
optical  systems.  This  correction  is  known  as  modulation  transfer 
function  compensation  (MTFC) ,  and  is  performed  by  the  convolution 
of  an  MTFC  matrix  with  the  scene  data.  Under  a  real-time  imple¬ 
mentation,  the  MTFC  could  be  modified  at  the  frame  rate.  This 
allows  real-time  control  of  image  sharpness  and  noise.  Another 
highly  desirable  feature  would  be  the  ability  to  modify  the  trans¬ 
formation  of  the  scene  intensities  under  the  real-time  control  of 
the  operator.  This  transformation,  known  as  dynamic  range  adjust¬ 
ment  (DRA) ,  can  be  implemented  by  the  use  of  a  look-up  table  to 
map  intensity  levels  in  any  appropriate  manner.  An  example  of 
its  use  would  be  to  bring  out  details  in  a  shadow  area. 


The  key  to  the  implementation  of  this  system  is  in  the 
approach  taken  to  handling  the  geometric  distortions.  The  solu¬ 
tion  to  the  problem  requires  both  the  use  of  a  computer  and  of 
special-purpose  real-time  hardware.  To  handle  a  variety  of 
taking  systems,  the  problem  must  be  cast  into  a  uniform  form  for 
use  by  the  special-purpose  hardware  by  using  software  within  the 
computer. 

The  first  step  is  to  generate  the  mapping  from  the 
object  surface  (r,s)  to  the  projection  surface  (a,£)  in  poly- 
nominal  form.  In  other  words,  Equations  C2-6)  and  (2-7)  must  be 
found  in  the  form: 

a  *  f'(r,s)  =  a1+a2r+a3s+a4r2+a5rs+a6s2+. . .  (3-1) 

£  *  f'(r,s)  =  b3+b2r+b3s+b4r2+b5rs+b6s2+. . .  (3-2) 

These  polynomial  coefficients  must  be  generated  by  the  computer 
software  and  be  passed  on  to  the  special-purpose  hardware  pre¬ 
vious  to  any  viewing  of  the  images.  In  this  manner  corrections 
for  any  type  of  detector  can  be  cast  into  a  single  form  which 
can  be  handled  by  special-purpose  hardware. 

In  order  to  determine  polynomial  coefficients,  a  grid 
of  points  must  be  mapped  from  the  object  surface  to  the  projec¬ 
tion  surface  using  the  appropriate  transformation  for  the  given 
detector  derived  from  the  initial  taking  parameters.  The  number 


of  grid  points  required  depends  upon  the  degree  of  the  polynomial. 
In  general,  the  number  of  points  used  must  be  greater  than  the 
number  of  coefficients  required  for  one  of  the  polynomials.  For 
example,  for  a  seventh  order  polynomial  requiring  36  coefficients, 
a  9x9  grid  of  81  points  should  be  sufficient.  Once  the  grid 
points  have  been  mapped  onto  the  projection  surface,  the  poly¬ 
nomial  coefficients  are  solved  for  by  the  least  squares  fitting 
which  minimizes  the  sum: 


1,3 


The  degree  of  the  polynomials  required  can  only  be  determined 
by  analysis  of  the  taking  parameters  for  the  systems  of  interest. 

These  geometric  correction  polynomials  permit  correc¬ 
tion  of  distortions  over  the  entire  extent  of  the  image.  Usually 
the  CRT  will  display  only  a  small  portion  of  this,  and  in  this 
local  area  it  is  possible  to  describe  the  distortions  by  a  much 
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smaller  polynomial.  It  is  suggested  that  for  viewing  purposes 
only  that  the  distortions  be  described  by  a  second  order  poly¬ 
nomial.  An  appropriate  CRT  resolution  for  displaying  local 
regions  of  the  image  might  be  512x512  pixels.  We  will  use 
this  resolution  for  example  purposes.  The  derivation  of  the 
appropriate  second  order  polynomial  from  the  higher  order  poly¬ 
nomial  must  be  evaluated  in  special-purpose  hardware.  It  is 
this  polynomial  evaluation  hardware  which  permits  the  transla¬ 
tion,  rotation,  and  magnification  of  the  image  by  the  operator 
at  frame  rates.  To  derive  the  second  order  polynomials,  a  change 
of  coordinates  is  performed.  We  explain  the  procedure  in  terms 
of  a  simple  translation  example.  It  is  desired  to  display  the 
sub-image  centered  at  the  point  (r0,sQ)  on  the  object  surface. 

A  transformation  of  coordinates  is  performed  which  maps  this 
point  into  the  point  (0,0).  This  coordinate  transformation  is 
given  by: 


r 


r-r. 


(3-4) 


The  correct  second  order  polynomials  are  obtained  by 
substitution  of  this  change  of  coordinates  into  the  polynomial 
correction  Equations  (3-1)  and  (3-2).  The  terms  out  to  second  order 
in  r'  and  s'  are  gathered  together  to  give  the  second  order  poly¬ 
nomials.  The  CRT  then  displays  the  interpolated  intensity  values 
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at  locations  corresponding  to  r'  going  from  -256  to  256  and  s" 
going  from  -256  to  256.  The  ability  to  rotate  and  magnify  or 
demagnify  the  sub-image  being  projected  is  obtained  in  the 
same  manner  by  the  use  of  the  appropriate  change  of  coordinate 
equations . 

After  the  second  order  polynomials  are  found,  the  next 
step  determines  where  each  of  the  512x512  points  to  be  displayed 
fall  in  projection  space  (the  space  in  which  the  original  measure 
ments  are  available) .  This  is  performed  by  a  polynomial  solver 
special-purpose  hardware  which  evaluates  the  second  order  poly¬ 
nomial  for  the  512x512  points  to  be  displayed.  Because  the  trans 
formation  is  now  only  of  second  order,  it  is  possible  to  design 
the  hardware  to  find  each  succeeding  point  by  using  only  addi¬ 
tions.  Reduction  of  the  polynomial  corrections  locally  to  second 
order  results  in  reducing  the  hardware  requirements  substantially 
at  this  point. 

The  final  special-purpose  hardware  required  is  a  two- 
dimensional  interpolator.  The  inputs  to  this  interpolator  in¬ 
clude  the  original  pixel  measurements  and  the  position  in  pro¬ 
jection  space  of  each  of  the  512x512  pixels  to  be  displayed. 

From  the  position  of  each  pixel  to  be  displayed  and  from  the 
original  intensity  measurements  in  the  immediate  vicinity,  the 
appropriate  intensity  values  to  be  displayed  are  derived. 
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Once  the  three  corresponding  points  have  been  measured 


by  the  operator,  the  images  can  now  be  brought  together  into 
registration.  This  requires  an  updating  of  the  distortion  poly¬ 
nomials.  These  corrections  must  be  calculated  in  the  computer. 

The  technique  uses  a  linear  transformation  of  both  the  images 
to  bring  them  together.  The  linear  transformation  is  of  the  form: 

z'  =  c^r+C2S+C2  (3-6) 

s'  =  c^r+Cj-S+Cg  (3-7) 


The  technique  used  to  derive  the  coefficients  finds  the  optimum 
amount  to  move  each  of  the  images,  given  the  errors  on  the  initial 
taking  parameters.  Updating  the  polynomial  corrections  requires 
a  substitution  or  change  of  variables  operation  similar  to  the 
technique  used  to  provide  for  the  translation  of  the  sub-image 
being  viewed. 


To  be  practical,  computer  aided  image  matching  must 
be  capable  of  processing  huge  amounts  of  image  data.  typical 
frame  might  come  from  an  original  photographic  image  which  is 
10X10  inches  square.  The  sampling  might  be  2000  points  per  inch 
The  image  data  for  this  case  is  an  array  of  20,000  by  20,000 
samples.  The  display  and  processor  must  be  capable  of  handling 
at  least  one  pair  of  these  images  and  should  be  capable  of  pro¬ 
viding  for  many  pairs. 

A  CRT  display  is  only  512x512  pixels  or  1024x1024. 
Larger  displays  could  be  developed,  perhaps,  but  are  not  really 
needed.  In  fact,  they  are  not  even  desirable.  If  we  had  a  CRT 
with  20,000  by  20,000  pixels  resolution,  we  would  need  a  micro¬ 
scope  for  most  image  examination  since  we  could  not  resolve  the 
image  detail  by  eye  alone.  Another  way  to  achieve  the  required 
magnification  on  such  a  super  CRT,  so  that  the  eye  can  perceive 
detail,  would  be  to  digitally  magnify  the  image  and  display  only 
a  small  image  portion  around  the  point  of  interest.  This  im¬ 
plies  that  many  CRT  pixels  are  used  to  display  each  original 
sample.  For  example,  for  10X  magnification,  the  full  screen 
would  be  used  to  display  an  area  only  2000  by  2000  samples.  One 
hundred  CRT  pixels  would  be  used  to  display  each  original  sample 
Of  course,  the  effect  of  magnification  by  digital  interpolation 
is  to  smooth  the  display  and  provide  optimum  resolution.  So 
each  original  sample  is  not  just  replicated,  but  on  the  average, 
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one  hundred  pixels  are  used  for  each  original  sample. 

This  is  just  like  examination  of  film  with  a  micro¬ 
scope.  The  eye  can  see  the  whole  frame,  but  cannot  resolve 
all  the  detail.  The  f ield-of-view  of  the  microscope  covers 
only  a  small  portion  of  the  film  image  area  and  provides  details. 

A  super  high  resolution  CRT  is  not  needed.  When  the 
entire  frame  is  displayed,  the  eye  can  not  resolve  fine  detail 
anyway.  A  512x512  or  1024x1024  CRT  can  provide  all  of  the  reso¬ 
lution,  detail  and  f ield-of-view  that  a  human  operator  can  in¬ 
terpret. 

To  make  a  normal  CRT  applicable  to  image  matching ,  the 
original  samples  must  be  processed  to  provide  displays  ranging 
from  full  frame  to  image  areas  small  enough  to  provide  full  reso¬ 
lution  and  detail.  Further,  this  processing  must  be  rapid  enough 
to  follow  the  operator's  manipulations.  Taking  clues  from  the 
use  of  light  tables  and  microscopes,  we  note  that  the  human  oper¬ 
ator  examines  the  whole  frame  or  large  portions  of  the  frame,  cen¬ 
ters  a  microscope  over  an  area  of  interest  and  observes  this  area 
in  greater  detail.  Often  the  operator  alternates  between  large 
and  small  scale  image  examination,  varying  magnification  and 
f ield-of-view  and  slewing  about  over  an  image.  Actually,  the 
f ield-of-view  of  a  zoom  microscope  varies  only  in  respect  to  the 
film  image.  The  virtual  image  presented  to  the  eye  is  of  fixed 
extent,  just  as  the  screen  of  a  CRT  is  of  fixed  size. 
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What  is  needed  for  practical,  computer  aided  manual 


image  registration  is  a  digital  processor,  which  can  handle  pairs 
of  20,000  by  20,000  sample  images  and  provide  arbitrary  slew 
and  zoom  for  a  CRT  by  full  image  alteration  during  the  l/30th 
second  framing  period.  Such  a  device  would  provide  all  of  the 
capability  of  a  light  table  and  zoom  microscope  coupled  with 
real-time  image  geometry  manipulation. 

This  digital  processor  and  display  for  image  matching 
is  well  within  the  current  state-of-the-art.  E-Systems  could 
build  this  image  manipulator  with  only  engineering  effort  -  no 
research  or  development  effort  would  be  needed. 
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4.0 


TEXTURE  MATCHING  CONCEPT  DEMONSTRATION 


In  this  section  the  concept  of  texture  matching  is 
demonstrated  with  the  use  of  digital  data  from  a  single  aerial 
image.  This  example  illustrates  texture  matching  as  an  extension 
of  conventional  correlation  matching  in  which  each  pixel  value  is 
first  replaced  by  some  measure  of  the  texture  in  its  immediate 
neighborhood . 

In  the  first  subsection,  a  sample  aerial  image  is  used 
to  investigate  a  conventional  correlation  technique  using  the 
normalized  cross  correlation  product.  The  object  here  is  to  derive 
a  second  image  which  no  longer  correlates  with  the  original  image. 
This  second  image  is  derived  by  modifying  the  original  image. 

Several  modifications  are  investigated  including  the  addition  of 
noise,  low  pass  filtering,  a  nonlinear  intensity  transformation, 
and  an  edge  enhancement  operation.  The  edge  enhancement  operation 
is  found  to  produce  an  image  which  no  longer  correlates  with  the 
original  image. 

In  the  next  subsection,  a  simple  texture  matching  approach 
of  finding  match  points  is  investigated.  The  original  image  and 
the  image  modified  by  the  edge  enhancement  operation  are  used. 

Each  pixel  in  both  images  is  replaced  by  a  measure  of  the  local 
activity  in  a  four  by  four  neighborhood  of  each  pixel.  After  this 
replacement,  the  normalized  cross  correlation  approach  is  found  to 
produce  reasonably  accurate  match  points.  Next,  factors  affecting 
the  accuracy  of  the  match  points  are  investigated.  Accuracy  is 
investigated  for  procedures  which  replace  every  other  pixel  in 
every  other  line  and  every  fourth  pixel  in  every  fourth  line  by 
the  same  measure  of  local  activity.  Using  these  alternative 

procedures,  roughly  the  same  accuracy  is  obtained. 
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Conventional  Correlation  Experiments 


Figure  4-1  shows  the  aerial  image  that  has  been  used 
in  Section  4  to  demonstrate  the  texture  matching  approach  to 
finding  image  match  points.  The  digital  representation  of  the  image 
has  been  obtained  by  scanning  the  original  picture  with  a  micro¬ 
densitometer.  The  digital  representation  of  the  image  consists  of 
approximately  1000  by  1000  pixels.  Each  pixel  is  represented  by 
a  ten  bit  reflectance  value  measured  by  a  densitometer.  The  digital 
image  has  been  reconstructed  to  form  the  picture  shown  in  Figure  4-1 
through  the  use  of  a  laser  scanning  device. 

In  this  subsection  we  would  like  to  find  modifications  of 
the  original  image  which  prevent  the  use  of  correlation  techniques 
to  find  match  points  between  the  original  and  the  modified  image. 
Hopefully,  these  modifications  will  be  representative  of  the  type  of 

problems  present  in  trying  to  match  dissimilar  images  of  the  type 

we  would  like  to  be  able  to  match.  Then,  in  the  next  subsection, 
the  texture  matching  extension  of  the  correlation  match  point 
technique  can  be  demonstrated  on  these  two  images. 

In  this  subsection  the  normalized  cross  correlation 
approach  is  used  to  determine  match  points  between  the  original 
image  and  the  modified  image  derived  from  the  original.  This 
is  a  powerful  approach  which  is  very  effective  for  determining 
match  points  between  similar  images.  For  convenience,  the  normalized 
cross  correlation  is  defined  here  for  the  one  dimensional  case. 

The  generalization  to  two  dimensions  is  straight  forward.  Let  x^ 
be  the  reference  window  which  we  would  like  to  match  with  some  area 

of  a  search  area  given  by  the  data  y^,  i.e. 
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FIGURE  4—1. 
SAMPLE  AERIAL  IMAGE 


reference  window: 
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This  normalized  cross  correlation  is  bounded  between  plus 

j 

and  minus  one.  High  correlation  corresponds  with  numbers  near  one,  j 

and  match  points  are  determined  by  searching  for  this  correlation 
peak.  If  it  is  necessary  to  search  large  areas,  one  problem  which 
can  arise  is  multiple  solutions.  In  general,  if  two  data  sets  differ 
in  their  intensity  distributions  by  only  a  multiplication  or  a  translation 
factor,  i.e. 

yi  =  Ax i  +  B  (4-6) 

then  their  cross  correlation  product  will  yield  exactly  one.  This 
can  be  demonstrated  by  substitution  into  the  above  defining  equations. 

This  is  an  important  property  which  permits  the  determination  of 
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match  points  for  a  much  broader  class  of  similar  images  than  would 
otherwise  be  true.  This  property  is  due  to  the  normalizing  proce¬ 
dure  used. 

4.1.1  No  Modification  Example 

Before  attempting  to  modify  the  sample  aerial  image, 
nominal  correlation  results  are  obtained  with  no  modification 
introduced.  The  general  procedure  followed  in  each  case  is  identical. 
A  reference  window  size  of  32  by  32  pixels  is  used.  A  search  area 
of  48  by  48  pixels  centered  over  the  reference  window  is  used. 

This  yields  a  17  by  17  element  correlation  space.  With  zero  error 
in  the  match  point  determination  process,  the  correlation  peak 
occurs  at  the  center  element  of  the  correlation  space  (the  ninth 
element  of  the  ninth  row) .  In  the  general  case,  the  reference 
window  is  taken  from  the  modified  version  of  the  aerial  image 
and  the  search  is  taken  from  the  unmodified  version  of  the  image. 

Ten  areas  have  '>een  chosen  at  random  from  the  image  for 
the  correlation  examples.  For  each  case  (i.e.  for  the  case  of 
no  modification  and  for  each  of  the  four  possible  modification  cases) 
the  identical  ten  areas  are  used.  In  other  words,  the  positions  of 
the  reference  windows  and  the  search  areas  are  unchanged. 

When  correlation  is  performed  without  the  introduction  of 
any  modification,  well  defined  correlation  peaks  are  obtained  for  each 
of  the  ten  examples.  Figure  4-2  shows  the  correlation  results  for 
two  of  the  ten  areas.  These  two  areas  are  referred  to  as  Area  A 
and  Area  B.  The  correlation  numbers  in  the  figure  have  been  scaled 
by  multiplying  by  a  thousand.  Thus  a  perfect  correlation  of  one 
appears  in  the  figure  as  1000.  For  ease  of  reading  the  center 
element  of  the  correlation  results  has  been  boxed  out  and  the 


correlation  peak  has  been  circled. 
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Area  A  and  Area  B  have  been  chosen  as  representative  of 
two  possible  extremes  present  in  the  data.  Area  A  represents  an 
area  which  appears  to  a  viewer  as  perfectly  flat,  devoid  of  any 
activity.  When  a  closer  look  is  taken  of  the  digital  data,  it 
becomes  apparent  that  a  noise  background  is  present  which  is  not 
apparent  to  the  viewer.  Because  of  this,  the  spectrum  of  frequencies 
present  in  the  reference  window  is  relatively  flat.  The  correlation 
peak  obtained  is  extremely  narrow.  This  is  due  to  the  relative 
strength  of  the  high  frequency  spectrum  components. 

Area  B  represents  an  area  of  the  picture  which  appears 
to  the  viewer  as  very  active.  Such  regions  are  found,  for  example, 
in  the  areas  in  which  the  cars  are  found.  The  correlation  peak  for 
this  area  is  much  broader  than  for  Area  A.  This  is  due  to  the  relative 
strength  of  lower  frequency  spectrum  components  present.  In  this 
case,  the  lower  frequency  components  are  much  stronger  than  the 
higher  frequency  components. 

4.1.2  Four  Modification  Examples 

We  would  like  to  find  modifications  of  the  original  aerial 
image  which  destroy  the  ability  of  conventional  correlation  techniques 
to  determine  match  points.  Toward  this  goal,  four  modifications  have 
been  tried. 

4. 1.2.1  Addition  of  Noise 

For  the  first  modification  example,  the  original  picture  was 
modified  with  the  addition  of  Gaussian  noise.  The  ten  reference 
windows  were  taken  from  the  modified  image  and  the  search  areas 
were  taken  from  the  original  image.  Figure  4-3  shows  the  results 
for  Area  A  and  Area  B  which  are  typical  of  the  results  for  all 
ten  cases.  The  standard  deviation  of  the  noise  introduced  was  +  40. 


AREA  A 


1 

2 

3 

4 

5 

6 

7 

9 

9 

10 

1 1 

_ \2_ 

1  3 

14 

15 

16 

17 

1 

7  j 

- .. 

-2l 

1 

24 

-7 

-20 

-29  - 

30 

-Sfc 

-8 

-2  3 

-40 

-12 

-29. 

27 

20 

z 

“  O  '• 

-  .  j 

-30 

57 

10 

0 

24 

-10 

SI 

-9 

-4  1 

.-64 

-24 

-7 

-44 

-20 

13 

3 

:  1 

13 

-22 

60 

49 

40 

18 

7 

-7 

6 

64 

a 

2 

-29 

-21 

-9 

4 

-  *3 

; \ 

- 12 

16 

-1 

26 

26 

-5 

1 1 

53 

57 

-  10 

-37 

-4 

22 

33 

-S 

5 

42 

-12 

-29 

0 

-1 

-1  6 

61 

-  19 

9  l 

45 

38 

77 

41 

95 

26 

6 

l  - 

-  25 

26 

-IS 

33 

71 

7  - 

13 

11 

3 

-26 

35 

33 

20 

30 

31 

7 

n 

74 

23 

54 

21 

38 

62 

55 

47 

35 

81 

82 

17 

31 

30 

41 

8 

•*>  t 

4  * 

14 

7 

43 

27 

S3 

105 

94 _ 

88 

71 

35 

-  19 

17 

37 

20 

63 

37 

9 

7 . 

69 

48 

17 

IS 

4 

117 

2  1  0 [8  47 j 214 

89 

-7 

-8 

14 

29 

79 

26 

10 

♦  4 

66 

3 

55 

47 

0 

46 

1  03 

78 

97 

83 

68 

17 

51 

SO 

33 

-14 

i  I 

36 

l 

34 

30 

91 

57 

S  2 

80 

so 

73 

40 

-  16 

SO 

38 

63 

90 

-5 

I  2 

20 

56 

-3 

4 

23 

-5 

-1  3 

10 

-5 

8 

57 

7 

-43 

17 

-15 

61 

25 

i  3 

14 

80 

7 

0 

11 

30 

89 

-9 

38 

-28 

9 

-29 

-29 

86 

25 

44 

-12 

1  4 

-4 

-8 

23 

-33 

-85 

-29 

6 

22 

-9 

0 

-6 

-35 

-IS 

1 

-32 

12 

0 

1  5 

IS 

-S 

-19 

-26 

-24 

-36 

-l  s 

-17 

3 

-34 

69 

52 

34 

-7 

-l 

-15 

-16 

J  b 

5 

1  2 

-  33 

-  25 

-46 

-30 

34 

-6 

19 

-I 

-13 

-ii 

-40 

i  1 

-12 

5 

-12 

1  7 

1 

9 

-26 

-  33 

-41 

-35 

43 

-24 

14 

-34 

-7 

16 

-14 

-53 

-56 

13 

48 

AREA  B 


1 

2 

3 

4 

5 

6 

7 

8 

9 

16 

1  1 

12 

13 

14 

IS  16  17 

n 

-  i  4 

33 

59 

93 

133. 

no 

240 

234 

218 

170 

112 

45 

-19 

-91  - 

133-142-155 

2 

-40 

22 

58 

98 

131 

161 

242 

2  65 

252 

218 

152 

75 

1 

-60- 

1 28-147-159 

.3 

-73 

4 

54 

99 

i  60 

2  10 

297 

359 

339 

296 

226 

1  24 

31 

-57- 

1 14-142-147 

4 

-77 

8 

72 

1  39 

2  1 1 

Z74 

35  3 

441 

4  15 

361 

276 

173 

79 

-10 

-77-108-1  10 

5 

-89 

1 

85 

1  6  ) 

2  42 

311 

406 

5  1  0 

4  94 

4  44 

359 

257 

160 

78 

7  -45  -52 

6 

-111 

-22 

83 

18  1 

2  7* 

362 

480 

593 

6  10 

5  49 

465 

359 

26  t 

1  66 

75  23  7 

7 

-  r  1 3 

*  1  2 

1  03 

217 

320 

431 

55  3 

671 

7  IS 

639 

54  8 

455 

354 

258 

169  93  48 

8 

-10  3 

1 

1  24 

245 

3S6 

469 

607 

729 

8  12 

7  SO 

670 

556 

433 

3o5 

178  81  22 

0 

-  b  3 

20 

1  40 

259 

380 

S  10 

67S 

823L980J  832 

699 

547 

405 

271 

ISO  50  -14 

1  0 

-39 

29 

1  40 

270 

394 

5  30 

660 

7  S  2 

8  15 

748 

626 

468 

346 

211 

93  -10  -75 

1 1 

2 

47 

14  2 

248 

336 

453 

S70 

673 

747 

706 

576 

441 

296 

165 

57  -28  -81 

1  2 

24 

34 

1  04 

213 

311 

4  22 

537 

621 

669 

653 

528 

399 

270 

151 

58  -16  -61 

13 

4  t 

42 

(0  2 

187 

265 

367 

469 

545 

S  89 

6  00 

489 

374 

266 

150 

61  -11  -57 

E 

41 

20 

72 

/  47 

226 

3  20 

41  6 

482 

516 

S  39 

44  9 

3  49 

238 

1  30 

46  -J9  -63 

36 

2S 

5  3 

1  31 

221 

294 

38  1 

429 

458 

465 

402 

299 

zos 

1  09 

3Z  -34  -34 

E 

25 

17 

35 

106 

201 

2S5 

31  8 

363 

3  8% 

3  87 

35  4 

2  49 

168 

77 

0  -64-1 lO 

17 

8 

32 

14 

75 

164 

2  11 

no 

30& 

3  33 

319 

296 

196 

1  J5 

40 

-39  -96-150 

FIGURE  4-3.  CORRELATION  EXAMPLE,  NOISE  MODIFICATION  ADDED 


43 


The  standard  deviations  of  the  reference  window  data  before  the 
addition  of  noise  were  +  63  and  +  198  for  Area  A  and  Area  E,  respectively. 

Inspection  of  the  two  correlation  peaks  shown  in  Figure 
4-3  shows  that  accurate  matches  can  still  be  obtained.  The  correlation 
peaks  maintain  approximately  the  same  width  and  are  correctly 
centered.  The  main  difference  appears  in  the  values  obtained  at 
the  maximum.  For  Area  A  and  Area  B  the  value  at  the  correlation 
peak  maximum  goes  from  1.0  to  0.847  and  0.980,  respectively.  The 
value  at  the  correlation  peak  maximum  is  dependent  upon  the  signal 
to  noise  ratio.  Area  B  represents  a  high  signal  to  noise  ratio 
case,  and  little  change  from  a  correlation  maximum  of  one  is 
observed.  Area  A  represents  a  lower  signal  to  noise  ratio  case, 
and  a  more  significant  change  from  one  is  observed  in  the  value  of 
the  correlation  peak  maximum. 

In  general,  the  standard  correlation  approach  to  matching 
appears  to  handle  the  addition  of  noise  very  adequately. 

4.1. 2. 2  Low  Pass  Filtering 

Next,  the  original  aerial  image  was  modified  by  the  use 
of  a  low  pass  filtering  operation.  A  simple  low  pass  filter  was 
used  which  makes  use  of  an  averaging  operation.  Each  pixel  value 
was  replaced  by  the  average  value  from  the  five  by  five  pixel  area 
immediately  surrounding  it.  For  each  of  the  ten  correlation  cases 
tried,  the  reference  windows  were  taken  from  the  modified  image  and 
the  search  areas  were  taken  from  the  original  image.  Figure  4-4 
shows  the  results  from  Area  A  and  Area  B  which  are  typical  of  the 
results  for  all  ten  cases. 
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In  all  cases  the  match  point  found  by  the  center  of  the 
correlation  peak  is  still  correctly  positioned.  However,  in  some 
cases  the  shape  of  the  peak  is  substantially  altered.  The  results 
for  Area  A  are  an  example  of  this.  With  the  addition  of  the  low  pass 
filtering,  the  peak  becomes  much  broader  and  has  a  maximum  at  a 
value  of  only  0.391.  Originally  in  this  case,  the  correlation  peak 
was  very  narrow  because  of  the  relative  strength  of  the  high  frequency 
spectrum  components.  The  low  pass  filtering  removes  the  high  frequency 
components  from  the  reference  window,  and  results  in  a  much  broader 
correlation  peak. 

For  Area  B,  the  original  correlation  peak  is  much  broader 
corresponding  with 'the  dominance  of  low  frequency  components.  The 
addition  of  the  filtering  operation  has  a  much  less  dramatic  effect 
on  the  correlation  peak.  It  becomes  slightly  broader,  and  peaks 
at  the  value  of  0.905. 

Again,  we  find  that  the  standard  correlation  approach  to 
matching  can  handle  modifications  introduced  through  the  use  of  a  low 
pass  filtering  operation. 

4. 1.2. 3  Nonlinear  Intensity  Transformation 

For  the  next  modification  tried,  a  nonlinear  transformation 
is  applied  to  each  intensity  value.  Figure  4-5  shows  the  nonlinear 
intensity  transformation  which  has  been  used.  With  this  transformation 
the  output  intensity  values  are  exponentially  related  to  the  input 
values.  For  each  of  the  ten  correlation  cases  tried,  the  reference 
windows  were  taken  from  the  modified  image  and  the  search  areas  were 
taken  from  the  original  image. 
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Again,  in  all  ten  cases  the  match  point  found  by  the 


center  of  the  correlation  peak  is  still  correctly  positioned.  In 
each  case  the  correlation  peak  becomes  slightly  broader  and  the 
value  at  the  maximum  is  slightly  less  than  one.  This  can  be  seen 
in  Figure  4-6  for  Area  1  and  Area  2.  In  these  two  cases  the 
correlation  peak  occurs  with  values  of  0.913  and  0.930,  respectively. 

Thus,  the  standard  correlation  approach  to  matching 
appears  to  handle  cases  of  reasonable  nonlinear  modifications 
applied  to  the  intensity  levels.  This  particular  property  is  due 
to  the  normalizing  procedure  used  in  the  cross  correlation  product 
we  have  used.  Of  course,  correlation  will  not  work  if  the  non¬ 
linear  intensity  transformation  is  dependent  upon  the  scene  content 
(the  situation  for  IR  imaging) . 

4 . 1 . 2 . 4  Edge  Enhancement  Operation 

The  final  modification  that  has  been  used  is  an  edge 
enhancement  operation.  Consider  the  three  dimensional  space  made 
up  of  the  two  spatial  dimensions  corresponding  with  the  pixel 
positions  and  the  intensity  values  corresponding  with  the  third 
dimension.  In  this  space  the  sequence  of  intensity  values  from 
the  original  image  define  a  surface.  For  each  pixel  position  the 
gradient  of  the  surface  gives  the  magnitude  and  direction  of  the 
maximum  slope  associated  with  each  pixel.  The  edge  enhancement  operator 
that  has  been  used  estimates  the  magnitude  of  gradient  associated 
with  each  pixel  position.  The  pixel  values  from  a  three  by  three 
neighborhood  surrounding  the  pixel  of  interest  are  used  in  making 
the  estimate.  For  the  pixel  position  (i,j)  let  the  neighborhood 
intensity  values  be  denoted  by: 
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FIGURE  4-6.  CORRELATION  EXAMPLE,  NONLINEAR  INTENSITY 
TRANSFORMATION  MODIFICATION  ADDED 
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Then  the  edge  enhancement  operator  that  has  been  used  is  defined  by: 


=  (x 
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i-l,j+l"Ai,  j+l+xi+l,j+l~*i-l,j-l'*i,  j-l"*i+l,  j-1 
sy  =  (xi-l, j-l+xi-l, j+xi-l, j+l“xi+l, j-l“xi+l, j"xi+l, j+1* /6  (4_8) 


-x. 


)/6  (4-7) 


Nsx2  +  s- 2 


i,j  "N~x  ~y 

s.  .  gives  an  estimate  of  the  magnitude  of  the  gradient  at  the  pixel 
1  •  J 

position  (i,j). 

There  are  two  ways  to  derive  this  estimate.  If  it  is 
assumed  that  the  nine  points  all  are  in  a  plane,  then  the  least 
squares  solution  for  the  plane  yields  the  above  solution. 
Interestingly  enough,  the  same  result  is  obtained  if  the  surface 
is  assumed  to  be  quadratic  in  its  dependence  on  position.  There 
are  many  edge  enhancement  operators  that  can  be  found  in  the  digital 
image  processing  literature.  The  above  estimate  gives  a  reasonable 
compromise  between  accurately  measuring  sharp  edges  and  providing 
some  resistance  to  the  effects  of  noise.  Using  fewer  pixel  values 
in  the  estimate  more  accurately  describes  sharp  slope  changes  but 
is  adversely  affected  by  the  presence  of  noise. 

For  each  of  the  ten  correlation  cases  tried,  the  reference 
windows  were  taken  from  the  modified  image  and  the  search  areas 
were  taken  from  the  original  image.  The  results  show  that  the 
application  of  the  edge  enhancement  operator  has  completely  destroyed 
the  effectiveness  of  the  correlation  approach  to  matching. 


(4-9) 
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Figure  4-7  shows  the  results  for  Area  A  and  Area  B.  Figure  4-8 
shows  the  distribution  of  the  position  of  the  correlation  maximum 
for  both  the  x  and  y  directions  for  the  ten  correlation  cases. 

These  distributions  are  approximately  uniform.  As  a  check  of  this 
basic  result,  the  larger  reference  window  size  was  used.  For  a  64  by 
64  pixel  reference  window  the  same  basic  result  was  obtained.  Match 
points  could  not  be  determined. 

This  result  is  exactly  what  has  been  sought  after.  Figure 
4-9  shows  the  edge  enhanced  image .  To  the  human  observer  it  is 
immediately  obvious  how  to  match  the  modified  image  to  the  unmodified 
one.  For  most  of  the  class  of  dissimilar  images  which  we  would  like 
to  match  this  is  the  case.  It  is  expected  the  differences  present 
between  these  two  images  that  prevent  the  correlation  approach  to 
matching  from  working  are  representative  of  the  differences  found 
in  many  of  the  dissimilar  images  that  we  would  like  to  match. 

4 . 2  The  Texture  Matching  Approach 

In  this  subsection,  a  simple  texture  matching  approach 
of  determining  match  points  is  demonstrated.  The  basic  idea  is  to 
develop  a  descriptor  that  would  be  associated  with  each  pixel  that 
would  be  a  measure  of  the  surround  of  that  pixel.  This  descriptor 
would  describe  the  local  scene  content  or  texture.  The  procedure 
used  is  to  replace  each  pixel  from  the  images  that  are  to  be 
matched  with  this  local  measure  of  texture.  In  general,  different 
replacement  algorithms  could  be  used  for  each  image.  Then  con¬ 
ventional  correlation  techniques  are  used  to  determine  match  points. 
In  the  example  presented  in  this  subsection,  the  normalized  cross 
correlation  product  is  used  for  the  correlation. 
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In  order  to  demonstrate  the  approach,  the  results  from 
the  previous  subsection  are  used.  Here,  it  was  demonstrated  that 
after  the  application  of  an  edge  enhancement  operation,  normal 
correlation  techniques  no  longer  are  successful  in  determining 
match  points.  This  is  true  despite  the  fact  that  a  human  observer 
would  have  no  trouble  in  matching  the  two  images.  Thus,  we  will 
attempt  to  match  the  original  aerial  image  with  the  image 
modified  by  an  edge  enhancement  operation  using  the  texture 
matching  approach. 

4.2.1  Texture  Matching  Example 

A  very  simple  quantity  is  used  here  as  a  measure  of 
local  texture.  For  each  pixel,  the  pixel  values  in  a  four  by  four 
neighborhood  surrounding  the  pixel  of  interest  are  used.  These 
pixel  values  are  averaged  to  find  a  mean  intensity,  and  then  the 
root-sum-square  deviation  about  the  mean  is  calculated.  Using 
the  same  notation  as  was  used  to  define  the  edge  enhancement 
operator,  this  transformation  is  defined  as: 
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This  measure  of  local  activity  or  texture  is  referred  to  as  the 
local  busyness  measure.  In  regions  near  edges  or  in  regions  of 
rough  texture  its  magnitude  is  large.  In  calm  or  flat  areas  its 
magnitude  is  small. 
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The  procedure  used  is  to  transform  both  the  original 
aerial  image  and  the  image  distorted  by  the  edge  enhancement 
operation.  Although  the  same  transformation  is  used  here  for  each 
image  of  the  pair,  in  general,  the  transformations  could  be  different. 
All  pixels  are  transformed  into  -the  edge  busyness.  Then  the  normalized 
cross  correlation  product  is  used  to  determine  match  points.  The 
reference  windows  are  selected  from  the  modified  image  and  the  search 
areas  are  selected  from  the  original  image.  Again  a  reference  window 
of  32  by  32  pixels  is  used  and  a  search  area  of  48  by  48  pixels  is  used. 
For  the  ten  matching  cases  that  are  tried,  exactly  the  same  areas  are 
selected  that  were  used  in  the  correlation  experiments  in  the  previous 
subsection. 

Figure  4-10  shows  the  correlation  results  for  Area  A  and 
Area  B.  As  can  be  seen,  a  well  defined  correlation  peak  appears  in 
both  cases.  The  center  of  the  peak  is  within  one  pixel  of  the  exactly 
correct  match  position.  Of  the  ten  examples  tried  this  is  true  in 
all  but  one  of  the  cases.  In  more  than  half  the  cases,  the  value  at 
the  correlation  peak  is  0.8  or  higher.  In  three  cases  the  value  at 

the  correlation  peak  is  approximately  0.55.  Thus,  the  transformation 

of  both  the  images  into  the  simple  texture  measure  has  permitted  the 
determination  of  match  points  in  most  cases  to  within  one  pixel  accuracy. 

Figure  4-11A  shows  the  texture  matching  results  for  the 
one  case  in  which  an  accurate  match  was  not  obtained.  In  this  case 

the  correlation  peak  consists  of  a  long  diagonal  ridge  with  values 

at  the  top  of  the  ridge  very  nearly  the  same.  In  this  direction 
the  position  of  the  match  is  poorly  determined  resulting  in  the  bad 
match  obtained.  To  obtain  a  match  on  this  region  of  the  image  it 
was  found  necessary  to  enlarge  the  reference  window  size.  Figure 
4-11B  shows  the  results  when  a  64  by  64  pixel  reference  area  is 

used.  A  correct  match  is  obtained. 
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A)  32  BY  32  PIXEL  REFERENCE  WINDOW 
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FIGURE  4-11.  TEXTURE  MATCHING  EXAMPLE,  FAILURE  CASE 
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4.2.2 


Resolution  Factors 


Several  areas  have  been  investigated  to  learn  about  some 
of  the  factors  affecting  the  resolution  obtained  in  using  the 
texture  matching  approach.  The  most  interesting  results  were 
obtained  when  different  pixel  replacement  resolutions  were  studied. 

In  the  previous  subsection  each  and  every  pixel  from  the  1000  by  1000 
pixel  images  were  replaced  by  their  local  busyness  measure  as 
defined  by  Equation  4-11.  This  replacement  resolution  is  referred 
to  as  one-by-one  replacement  resolution.  Alternately,  it  is  possible 
to  replace  only  every  other  pixel  in  every  other  line.  This  gives 
a  500  by  500  local  busyness  representation  of  the  original  images. 

This  is  referred  to  as  two-by-two  replacement  resolution.  Replacing 
every  fourth  pixel  in  every  fourth  line  gives  four-by-four  replacement 
resolution. 

It  is  important  to  distinguish  what  is  being  done  here 
from  a  similar  but  different  process.  One  could  take  the  original 
image  and  abstract  every  other  pixel  from  every  other  line.  Then 
the  local  busyness  could  be  determined  using  this  subset  of  the 
original  pixels.  This  is  not  what  is  being  done  here.  Rather, 
local  busyness  is  just  being  calculated  for  every  other  pixel 
in  every  other  line  exactly  as  orginally  defined. 

The  use  of  two-by-two  and  four-by-four  replacement 
resolution  have  been  studied.  Again  ten  matches  are  attempted 
for  both  resolutions.  In  each  case  the  same  areas  of  the  images 
are  used  for  the  reference  windows  and  for  the  search  areas.  For 
two-by-two  replacement  resolution,  the  original  32  by  32  pixel 
reference  window  becomes  16  by  16  elements  wide,  the  original  48 
by  48  pixel  search  area  becomes  24  by  24  elements  wide,  and  the 
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size  of  the  correlation  space  becomes  9  by  9  instead  of  17  by  17 
For  the  four-by-four  replacement  resolution,  the  original  32  by 
32  pixel  reference  windows  becomes  8  by  8  elements  wide,  the 
original  48  by  48  pixel  search  area  becomes  12  by  12  elements  wide, 
and  the  size  of  the  correlation  space  becomes  5  by  5. 

Figure  4-12  shows  the  matching  results  for  two-by-two 
replacement  resolution  for  Area  A  and  Area  B.  For  a  correct  match 
the  correlation  peak  should  occur  at  the  fifth  column  of  the  fifth 
row.  Each  row  or  column  off  this  corresponds  with  a  two  pixel 
displacement.  Figure  4-13  shows  the  matching  results  for  four-by- 
four  replacement  resolution  for  Area  A  and  Area  B.  For  a  correct 
match  the  correlation  peak  should  occur  at  the  third  column  of  the 
third  row.  Each  row  or  column  off  this  corresponds  with  a  four 
pixel  displacement. 

We  would  like  to  compare  the  resolution  in  match  point 
determination  for  the  three  replacement  resolution  approaches 
(one-by-one,  two-by-two,  and  four-by-four) .  In  order  to  obtain 
the  position  of  the  correlation  peak  to  the  accuracy  of  a  fraction 
of  a  pixel,  a  simple  fitting  procedure  has  been  used.  The  x  and 
y  directions  are  treated  independently.  .  In  each  direction  a  three 
point  fit  to  a  parabola  is  performed  to  determine  the  peak  position. 

Table  4-1  summarizes  the  results  of  this  procedure.  In 
most  cases,  the  accuracy  of  the  match  is  correct  to  within  a  fraction 
of  a  pixel  in  both  directions.  For  one  case  (case  5) ,  a  reasonable 
match  was  not  obtained  and  this  case  has  not  been  used.  For  the 
nine  cases  left  the  average  accuracy  in  the  x  and  y  directions 
has  been  obtained  from  the  square  root  of  the  sum  of  the  errors 
squared.  The  average  radial  error  in  the  match  point  determination 


AREA  A 


1 

2 

3 

S 

6 

7 

a 

_ 2. 

1 

-I2S- 

100 

-SS 

-96- 

•104 

-6  2 

10 

73 

89 

z 

-  1 

56 

9  { 

53 

2  S 

57 

1  (9 

LOS 

45 

3 

16 

(42 

209 

236 

274 

23S 

170 

l  1  t 

S4 

4 

-60 

1  22 

295 

467 

613 

45  2 

178 

64 

39 

5 

-20  2 

-30 

26 1 

578 

pssj 

S3  1 

ISO 

1 

-LB 

6 

-239- 

t  19 

126 

390 

510 

323 

20- 

113 

-79 

1 

-110- 

102 

1  7 

90 

6S 

0 

-87- 

1  46 

-71 

9 

49 

-37 

-36 

-so- 

131- 

123 

-  33 

30 

91 

4 

135 

-18 

-96- 

136- 

199- 

147 

-14 

91 

167 

AREA  B 


1 

2 

3 

4 

5 

6 

7 

8 

9 

1 

345 

441 

4S7 

S38 

5  35 

439 

411 

413 

249 

2 

368 

479 

49  1 

608 

625 

529 

473 

492 

328 

3 

40  S 

516 

515 

681 

73  2 

614 

493 

520 

396 

4 

441 

528 

50  l 

702 

816 

7  02 

530 

552 

443 

5 

442 

SOI 

479 

7 1  a  | 

1 744 

47  1 

461 

376 

6 

395 

467 

415 

642 

843 

699 

39  3 

352 

305 

7 

34  9 

437 

37  6 

541 

753 

65  2 

372 

269 

22  4 

8 

325 

426 

373 

472 

695 

65  3 

41  1 

25  6 

195 

9 

273 

409 

373 

421 

632 

64  0 

451 

270 

195 

Figure  4-12.  TEXTURE  MATCHING  EXAMPLE,  TWO-BY-TWO  PIXEL 
REPLACEMENT  RESOLUTION 
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TABLE  4-1 


SUMMARY,  MATCH  POINT  ERRORS,  TEXTURE  MATCHING  EXAMPLES 


^^'''X^RR  3) 
Case 

Ax 

Ay 

4X4 

2X2 

1X1 

4X4 

2X2 

1X1 

1  (A) 
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-0.12 

0.02 
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0.27 

0.00 
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-0.84 

7 

-0.67 
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a)  PRR  =  Pixel  Replacement  Resolution 
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is  then  obtained  by  adding  in  quadrature  the  average  error  in  the 
x  and  y  directions.  The  results  yield  +  1.04,  +  1.04,  and 
+  1.08  for  the  4  x  4,  2  x  2,  and  lxl  replacement  resolutions, 
respectively.  In  other  words,  all  three  replacement  resolutions 
determine  the  match  point  to  approximately  the  same  accuracy. 

This  is  an  important  result.  The  use  of  the  four-by-four 
replacement  resolution  requires  one-sixteenth  fewer  calculations 
than  the  one-by-one  replacement  resolution  case.  This  result  is 
not  completely  unexpected.  The  calculation  of  local  busyness 
from  one  pixel  to  the  next  is  highly  redundant.  However,  it  was 
expected  that  there  would  be  some  fall  off  in  match  point  resolution, 
in  going  from  two-by-two  to  four-by-four  replacement  resolution. 

One  other  factor  affecting  match  point  resolution  was 
investigated.  The  effect  of  increasing  the  reference  window  size 
from  32  by  32  to  64  by  64  pixels  was  studied.  It  was  found  that 
the  match  point  resolution  was  improved  by  approximately  a  factor 
of  two. 

4 . 3  Some  General  Conclusions 

A  few  general  conclusions  can  be  drawn  from  these  simple 
experiments  conducted  to  demonstrate  the  texture  matching  approach 
to  matching  dissimilar  images.  In  most  parts  of  most  images,  the 
power  present  in  the  lower  frequency  components  far  exceeds  that 
found  in  the  higher  frequency  components.  Using  the  conventional 
approach  to  correlation  matching,  this  causes  the  low  frequency 
components  to  dominate  the  selection  of  the  match  point.  In 
dissimilar  images  of  the  type  we  would  like  to  match,  it  is  usually 
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found  that  the  lower  frequency  components  from  one  image  to  the 
next  do  not  agree.  Thus,  to  have  any  hope  of  success  in  matching 
dissimilar  images,  the  technique  must  not  use  information  related 
to  the  low  frequency  behavior  of  the  images .  Note  that  in  using 
the  texture  matching  approach  that  has  been  suggested,  this  will 
always  be  the  case.  If,  in  general,  each  pixel  is  replaced  by 
some  reasonable  measure  of  local  texture,  then  the  low  frequency 
components  of  the  image  are  effectively  discarded. 

Next  consider  the  size  of  the  neighborhood  used  to 
determine  the  local  texture  measure.  In  the  previous  subsection, 
a  four-by-four  pixel  neighborhood  was  used  to  determine  the  local 
busyness  texture  measure.  What  would  be  the  effect  on  match  point 
resolution  if  the  size  of  this  neighborhood  is  made  larger?  In 
the  last  subsection  we  demonstrated  that  it  is  not  necessary  to 
replace  every  pixel  by  the  texture  measure.  Further,  not  replacing 
every  pixel  does  not  effect  the  match  point  resolution  obtained. 

This  is  due  to  the  high  correlation  present  in  the  calculations  of 
the  texture  measures  from  one  pixel  to  the  next.  As  the  size  of  the 
neighborhood  used  is  increased  this  correlation  should  become  even 
stronger.  Thus,  it  is  suggested  that  as  size  of  the  neighborhood 
used  is  increased,  the  match  point  resolution  obtained  will  corre¬ 
spondingly  increase  (get  worse) .  For  this  reason,  it  is  suggested 
that  the  size  of  the  neighborhood  used  in  determining  the  texture 
measure  be  limited  to  four-by-four ,  eight-by-eight,  or,  perhaps, 
sixteen-by-sixteen  pixels  about  the  pixel  of  interest.  In  some 
cases  larger  areas  may  be  needed,  but,  in  general,  it  is  expected 
that  the  optimum  area  size  will  be  quite  small.  It  is  further 
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suggested  that  as  a  general  rule  that  the  pixel  replacement 
resolution  to  be  used  correspond  with  half  the  distance  used  in 
defining  the  neighborhood.  In  other  words,  for  a  four-by-four 
neighborhood,  a  two-by-two  pixel  replacement  resolution  be  used 
or  for  a  eight-by-eight  neighborhood,  a  four-by-four  pixel 
replacement  resolution  be  used.  Closer  spaced  pixel  replacement 
resolution  would  only  increase  the  computational  load  with  no 
increase  in  match  point  resolution.  Broader  spaced  pixel  replace¬ 
ment  resolution  must  eventually  degrade  the  match  point  resolution 
Time  permitted  for  this  program  did  not  permit  a  detailed  study 
to  confirm  these  suggestions. 
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In  this  section  refinements  to  the  basic  texture 
matching  approach  to  image  matching  are  discussed.  Algorithms 
are  suggested  that  could  be  coded  in  FORTRAN  IV  and  tested  on 
the  DIAL  facility  at  USAETL.  Before  discussing  these  refinements, 
the  basic  approach  is  summarized. 

The  two  images  to  be  matched  are  first  corrected  for 
all  known  geometric  distortions  using  the  initial  knowledge  of  the 
taking  conditions  (see  Section  2) .  As  part  of  this  process,  re¬ 
sampling  of  the  images  is  performed  if  necessary.  This  resampling 
process  aligns  the  coordinate  axes  used  in  the  two  images  and 
provides  for  equal  sample  spacing  for  the  two  images.  It  is  assumed 
that  the  dominant  effect  of  inaccuracies  in  the  taking  conditions 
give  a  translational  mismatch  between  the  two  images.  Each  of  the  images 
is  independently  transformed  on  a  pixel  by  pixel  basis.  Each  pixel  intensity 
value  is  replaced  by  a  local  measure  of  the  local  scene  content  (local  S 
texture) .  A  reference  window  is  chosen  in  one  image  and  a  larger 
search  area  is  chosen  in  the  other  image.  The  match  point  is  j 

determined  using  the  normalized  cross  correlation  product  (see 

( 

Section  4.1).  The  match  point  is  determined  by  the  position  of  the 

i 

correlation  peak.  A  threshold  can  be  set  using  the  value  of  the 
correlation  parameter  at  the  peak  to  determine  if  a  good  match 
has  been  found.  This  process  can  be  repeated  across  the  two  images 
to  obtain  match  points  as  frequently  as  desired. 

In  the  first  subsection  a  few  simple  transformations 


are  discussed  for  image  matching.  In  the  next  subsection  texture 
measures  which  can  be  derived  from  local  spectral  analysis  are 
discussed.  The  local  spectral  analysis  is  performed  with  the  use 


of  the  two  dimensional  discrete  Fourier  transform.  This  is  a 
powerful  approach  and  generates  a  broad  class  of  texture  measures. 

It  becomes  clear  in  this  subsection  that  to  describe  texture  in 
detail  requires  more  than  just  a  single  texture  measure.  In  the 
next  subsection  the  texture  matching  technique  is  formally  defined, 
and  this  technique  is  generalized  to  handle  multiple  measures  of 
texture.  The  next  subsection  discusses  a  variety  of  special 
problems  which  can  interfere  with  texture  matching.  Methods  are 
discussed  for  dealing  with  data  noise  characteristics,  different 
detector  spatial  acceptances,  drastically  different  intensity 
distributions,  and  coherency  effects  found  in  dealing  with  imaging 
radars. 

5.1  Some  Simple  Transformations 

The  two  transformations  discussed  in  this  subsection  are 
not  strictly  speaking  transformations  which  measure  local  texture. 
However,  they  are  transformations  which  allow  image  matching  between 
a  broader  class  of  images  than  would  be  true  for  conventional 
correlation  matching.  The  two  transformations  replace  the  original 
image  by  either  a  high  pass  frequency  filtered  version  of  the 
image  or  an  edge  enhanced  version  of  the  image. 

5.1.1  High  Pass  Filtering 

In  Section  4.3  it  was  pointed  out  that  in  most  of  the 
dissimilar  images  of  the  type  we  would  like  to  match,  the  lower 
frequency  components  do  not  agree.  Therefore,  one  approach  to 
the  problem  is  to  remove  these  components  through  the  use  of  two 
dimensional  high  pass  filtering.  An  example  of  where  this  approach 
might  be  helpful  would  be  two  overhead  images  taken  from  the  same 
perspective  but  at  different  times  of  the  day.  An  appropriate  high- 


pass  filter  can  be  developed  with  standard  signal  analysis  tools. 
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5.1.2 


Edge  Enhancement  Operation 

A  slightly  broader  class  of  images  could  be  matched  if 

images  are  first  processed  by  an  edge  enhancement  operator  such  as 

was  defined  in  Section  4. 1.2. 4.  This  operator  measures  the  magnitude 

of  the  gradient  of  the  intensity  level  surface.  Notice  that  for  an 

operator  of  this  type,  an  edge  that  goes  from  a  low  intensity  level 

to  a  high  intensity  level  will  yield  the  same  response  as  an  edge 

that  goes  from  a  high  intensity  level  to  a  low  intensity  level. 

This  is  because  the  magnitude  of  the  gradient  is  used.  In  this  way 

it  differs  from  the  use  of  a  high  pass  filter  and  permits  a  slightly 

broader  class  of  images  to  be  matched.  For  actual  implementation  it 

is  suggested  that  the  specific  edge  enhancement  operator  defined  in 

Section  4. 1.2. 4  be  used.  As  was  pointed  out,  this  operator  provides 

a  reasonable  compromise  between  accurately  measuring  sharp  edges 

and  providing  some  resistance  to  the  effects  of  noise. 

5 . 2  Texture  Measures  from  Spectral  Analysis 

In  this  subsection  texture  measures  which  can  be  derived 

from  local  spectral  analysis  are  considered.  For  each  pixel,  a 

neighborhood  is  defined  about  which  you  wish  to  measure  parameters 

associated  with  the  local  texture.  For  the  neighborhood,  the  local 

spectral  components  are  determined  using  the  two  dimensional  discrete 

Fourier  transform.  For  computational  simplicity,  the  neighborhood 

size  dimensions  are  usually  chosen  to  be  some  power  of  two  pixels 

wide.  For  example,  the  neighborhood  size  might  be  four-by-four, 

eight-by-eight,  or  sixteen-by-sixteen  pixels  wide.  Once  the 

frequency  components  are  determined,  the  texture  measure  is 

determined  from  these  component  values.  The  general  approach 

here  uses  the  square  root  of  the  power  present  in  similar  groups  of 

frequency  components.  For  example,  the  power  present  in  the  low  frequency 
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terms  or  the  high  frequency  terms  could  be  used.  It  quickly  becomes 
evident  that  to  describe  texture  in  detail  requires  more  than  one 
texture  measure.  In  order  to  accommodate  more  than  one  texture 
parameter,  the  general  matching  procedure  that  has  been  presented 
needs  to  be  generalized.  This  is  discussed  in  Section  4.3.2. 

To  be  able  to  group  together  similar  spectral  components 
in  a  reasonable  manner  requires  some  knowledge  of  the  properties 
of  the  discrete  Fourier  transform.  These  properties  are  discussed 
here  first  using  the  one  dimension  discrete  Fourier  transform. 

Given  a  sequence  of  N  equally  spaced  intensity  values 
denoted  by  xQ,  x,,  ...»  xN_i»  its  discrete  Fourier  transform 
is  given  by 


1  N-l  -2Tri  ju 

X  =  —  Z  x .  e  N 

u  «sTf  j=0  3 

where  i  =  iTT.  The  inverse  transform  is  given  by 

,  N-l 


(5-1) 


Xj  = 


2iri  ju 
Xu  e  N 


(5-2) 


\|  N  U=0 

It  should  be  noted  that  there  is  not  a  universally  accepted  form  of 
these  equations.  Some  times  a  reversal  of  signs  is  used  in  the 
kernels.  Some  authors  prefer  to  place  all  scaling  constants  in 
the  inverse  transform  equation.  The  normalization  used  here  leads 
to  a  particularly  simple  form  for  Parseval's  relationship. 


N-l 

Z 

xi  * 

N-l 

Z 

i=0 

u=0 

(5-3) 
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Note  that  the  spectral  components  as  defined  are  complex. 

At  first  glance  it  appears  that  from  N  independent  spatial  components 
we  have  produced  2N  independent  frequency  components.  Actually  this 
is  not  true.  Some  of  the  components  are  only  real.  Others  are 
closely  related.  The  relationship  between  spectrum  components  is 
most  easily  expressed  if  the  periodic  nature  of  the  spectrum 
components  is  understood.  Normally,  the  frequency  components  for 
u  =  0,  1,  . ..,  N-l  are  calculated.  This  is  not  the  only  choice. 
Equation  5-1  has  the  property  that  the  output  is  periodic  in  N,  i.e. 


X  =  x  „ 
u  u+N 


(5-4) 


for  any  u.  It  is  convenient  for  our  purposes  to  calculate  the 
frequency  components  for  u  =  -N/2+1,  . ..,  -1,  0,  1,  ...»  N/2.  In 
terms  of  these  components  the  relationships  that  reduce  the  number 
of  independent  frequency  components  to  N  are: 

Im(XQ)  =  Im(XN/2)  =  0  (5-5) 

Xu  =  X_u  ,  u  =  1,  2,  .  ..,  N/2-1  (*  =  conjugate)  (5-6) 


This  particular  choice  of  the  range  of  u  used  in  finding 
the  frequency  components  is  very  important  to  the  discussion  to 
follow.  With  this  choice,  the  frequency  of  the  wave  associated  with 
the  frequency  component  Xu  is  proportional  to  the  absolute  value  of  u. 
The  spectral  components  X2  and  X_2  are  associated  with  a  wave  of  twice 
the  frequency  of  components  X^  and  X_^,  etc.  We  will  be  using  the 
two  dimensional  generalization  of  the  discrete  Fourier  transform 
to  produce  frequency  components  Xy  With  this  particular  choice 
of  the  range  on  u  and  v,  it  is  relatively  easy  to  identify  closely 
related  frequency  terms.  With  the  more  conventional  choiee  this 
task  becomes  very  awkward. 


The  use  of  the  discrete  Fourier  transform  can  be  extended 

in  a  straigth  forward  manner  to  two  dimensions.  In  practice  the 

two  dimensional  form  is  calculated  by  a  sequence  of  applications 

of  the  one  dimensional  transform.  The  input  data  is  arranged  in 

a  square  matrix.  Each  column  is  transformed  using  the  one  dimensional 

transform.  The  rows  of  the  resulting  matrix  are  then  each  transformed 

using  the  one  dimensional  transform.  The  output  are  the  two  dimensional 

frequency  components  X„  .  Alternately,  the  rows  can  be  transformed 

first,  followed  by  the  columns. 

The  transform  can  be  implemented  using  the  well  known 

fast  algorithm,  the  fast  Fourier  transform  (FFT) .  Using  the  FFT, 

the  number  of  additions  and  multiplications  become  proportional  to 

2 

NlogN  instead  of  N  .  For  two  dimensions,  the  number  of  additions 

2  4 

and  multiplications  is  proportional  to  2N  logN  instead  of  N  * 

The  general  prescription  used  to  develop  texture  measures, 
groups  together  like  frequency  components.  These  are  combined 
to  give  a  texture  measure  using  the  general  prescription 


where  the  sum  over  u,v  includes  all  the  like  frequency  components 
that  have  been  selected.  Note  that  with  this  prescription  all 
phase  information  is  discarded. 

It  is  interesting  to  point  out  that  the  local  busyness 
texture  measure  is  of  this  form.  This  is  the  texture  measure 
that  was  used  in  Section  4.2  to  demonstrate  the  texture  matching 
approach.  This  measure  is  in  fact  equal  to  the  square-root  of  AC 
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energy  contained  within  the  block.  In  other  words,  to  construct 
this  texture  measure  (tB) ,  the  above  sum  is  taken  over  all  u,v 
except  (u , v )  =  (0,0).  The  equivalence  of  this  definition  of 
local  busyness  to  the  definition  used  in  Section  4.2.1  can  be 
demonstrated  using  Parseval's  relationship  (Equation  5-3). 

It  should  be  noted  that  the  term  Xq  q  is  never  used 
in  defining  any  of  the  texture  measures.  This  is  the  term  which 
is  proportional  to  the  average  intensity  value  of  the  block. 

Texture  measures  defined  in  this  manner  can  be  combined 
to  form  new  measures.  For  example,  suppose  the  texture  measures 
t^  and  ty  have  been  defined.  Here  t^  is  a  measure  of  strength 
of  frequency  components  which  are  only  a  function  of  the  x  direction, 
ty  is  a  measure  of  the  strength  of  frequency  components  which  are 
only  a  function  of  the  y  direction.  Then  a  new  texture  measure 
tx/y  could  be  defined  by  dividing  tx  by  ty.  This  would  give  you  a 
relative  measure  of  the  strength  of  waves  which  are  a  function  of 
the  x  direction  only  to  waves  which  are  a  function  of  the  y  direction 
only.  Alternately,  tx  and  ty  could  be  normalized  by  the  local 
busyness  measure,  i.e. 

fcX  =  tX//fcB'  fcY  =  tY//fcB  (5-8) 

In  the  X  direction,  for  example,  this  would  give  a  ratio  of  the 
strength  of  waves  which  are  a  function  of  the  x  direction  only 
to  the  overall  strength  of  all  waves  present. 

First,  specific  texture  measures  are  suggested  for 
neighborhoods  defined  to  be  four-by-four  or  eight-by-eight  pixels 
wide.  Then,  a  refinement  is  suggested  for  calculating  the  frequency 
components.  The  procedure  suggested  helps  to  reduce  unwanted  effects 
generated  by  the  very  strong  low  frequency  components  present  in 
most  images. 
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5.2.1  Four-by-Four  Pixel  Neighborhood 

For  a  four-by-four  neighborhood,  the  sixteen  complex 
frequency  coefficients  are: 

X2,-l  X2,0  X2,l  X2 , 2 

Xl,-1  X1 , 0  Xl,l  X1 , 2 

X0,-l  X0 , 0  X0,l  X0 , 2 

X-l,-l  X-1 , 0  X-l,l  X-1 , 2 

In  this  form  it  is  easy  to  group  together  similar  frequency  terms 
and  form  texture  measures. 


For  a  first  example,  the  low  or  high  frequency  terms  can 


be  grouped  together  to  form  a  measure  of  the  strength  of  the  low 
or  high  frequency  terms.  These  are  defined  by 


The  procedure  used  is  a  simple  case  of  a  general  procedure  of 
grouping  together  frequency  components  using  rings  defined  about 
(0,0)  in  frequency  space.  This  is  further  illustrated  in  the 
next  subsection. 


Alternatively,  the  strength  of  waves  which  are  only  a 
function  of  the  x  direction,  the  strength  of  waves  which  are 
only  a  function  of  the  y  direction,  and  the  strength  of  waves  which 
are  a  function  of  both  directions  can  be  measured.  These  are 
defined  by: 


76 


k-ij 

l2+M 
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^  (5-11) 

«Y  4" 

fcXY  =N> 
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xi„ 

J2  (5-12) 

x-i, 

-4 

»-] 

2+  : 

(5-13) 


5.2.2  Eight-by-Eight  Pixel  Neighborhood 

For  an  eight-by-eight  neighborhood,  there  are  64  complex 
frequency  coefficients.  These  are  symbolically  represented  in  Figure 
5-1.  In  this  figure  the  low,  medium  and  high  frequency  terms 
have  been  grouped  together  using  rings  superimposed  on  the  frequency 
components.  The  rings  are  defined  by'vTu2  +  v2  =  0.5,  1.5,  and  2.5. 
The  components  used  to  define  the  strength  of  the  low,  medium  and 
high  frequency  components  are  listed  in  Table  5-1.  These  three 
texture  measures  are  denoted  by  tL,  tM,  and  tg. 

Figure  5-2  shows  alternative  groupings  of  similar 
frequency  terms.  Here  pie  shaped  regions  about  the  (0,0)  element 
are  used  to  define  three  texture  measures.  These  are  referred  to 
as  tx,  ty»  and  tXy.  The  components  used  in  each  of  these  texture 
measures  are  listed  in  Table  5-1. 

These  six  measures  should  be  more  than  sufficient  to 
describe  the  texture  of  an  eight-by-eight  neighborhood.  Endless 
alternative  descriptions  could  be  formed  by  using  ratios  of  the 
measures  or  by  normalizing  by  the  total  AC  energy  (tg) .  In  any 
case,  it  should  not  be  necessary  to  use  more  than  a  half  dozen  or 
so  measures. 
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LOW  FREQUENCY  REGION 


MEOIUM  FREQUENCY  REGION 

HIGH  FREQUENCY  REGION 


AS090480-2I1 


FIGURE  5-1.  DIVISION  OF  FREQUENCY  COMPONENTS  DEFINING 

THE  LOW,  MEDIUM,  AND  HIGH  FREQUENCY  TEXTURE 
MEASURES 


TABLE  5-1 
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If  it  was  desired  to  describe  texture  more  finely,  the  frequency 
components  could  be  simultaneously  divided  by  the  radial  and  pie 
shaped  regions  shown  to  give  nine  texture  components . 

Note  that  to  define  texture  measures  by  grouping 
together  similar  frequency  terms,  it  has  not  been  necessary  to 
know  the  N  relationships  which  reduce  the  64  complex  terms  to  64 
independent  numbers.  These  relationships  are  useful,  however,  for 
the  implementation  of  the  calculation  of  the  frequency  components. 
For  convenience  in  expressing  these  relationships,  N/2  is  denoted 
by  M.  These  relationships  are: 


Im(x0i„)  =  In.(xM/0)  = 

Im(X0.M>  '  *  0 

(5 

* 

X  =  X  , 

u,v  =  0,1,  . . .M-l 

(5 

u,v  -u,-v 

* 

* 

X„  M  =  X  „  M  and  XM  „ 
U,M  -U,M  M,  U 

=  XM,-u  u  =  1/2,  ...,M-1 

(5 

The  use  of  these  relationships  reduces  the  number  of  one  dimensional 
discrete  Fourier  transforms  that  must  be  performed  from  2N  to 
approximately  3/2  N.  Figure  5-3  helps  illustrate  this  for  the 
eight-by-eight  neighborhood  case.  This  figure  shows  an  example 
of  a  minimum  set  of  components  that  must  be  explicitly  calculated. 
The  remaining  components  can  be  found  from  the  above  relationships. 
Eight  one  dimensional  transforms  would  be  required  on  the  eight 
columns,  but  only  five  one  dimensional  transforms  would  be  required 


on  the  rows  instead  of  eight. 


<•)  REAL  FREQUENCY  COMPONENT 

•  COMPLEX  FREQUENCY  COMPONENT 
Miiil  'NOEPENDENT  SET  OF  FREQUENCY  COMPONENTS 
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FIGURE  5-3. 


DEFINITION  OF  AN  INDEPENDENT  SET 
OF  FREQUENCY  COMPONENTS 


5.2.3 
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Low  Frequency  Background  Subtraction 
In  most  images  the  very  lowest  frequency  components 
have  by  far  the  most  power.  Between  dissimilar  images  of  the 
type  to  be  matched,  the  low  frequency  behavior  from  one  image 
to  the  next  is  usually  very  different.  Because  of  their  relative 
strength  and  because  of  the  properties  of  the  discrete  Fourier 
transform,  these  very  low  frequency  components  can  strongly  interfere 
with  the  texture  measures  which  have  been  defined  in  this  subsection. 
This  problem  is  illustrated  in  Figure  5-4.  Here  we  consider  what 
happens  when  a  ramp  signal  is  texture  analyzed.  From  a  texture 
point  of  view  the  signal  is  essentially  flat,  devoid  of  any  texture. 
Thus,  the  texture  measures  should  give  null  readings. 

Unfortunately,  this  is  not  what  happens.  Assuming  a 
block  size  of  four  pixels,  the  signal  is  divided  into  blocks  and 
each  block  is  analyzed  separately.  The  Fourier  transform  of  the 
block  gives  the  spectrum  of  the  discrete  periodic  signal  generated  by 
infinitely  replicating  the  four  pixel  values  in  the  block.  The  periodic 
nature  of  the  discrete  sequence  represented  by  the  discrete  Fourier 
transform  can  be  shown  from  Equation  (5-2) .  The  discrete  periodic 
sequence  resulting  is  shown  in  Figure  5-4.  Calculation  of  the 
spectral  components  gives  non- zero  high  frequency  terms,  and  as  a 
result  the  texture  of  the  block  is  misrepresented. 

The  solution  of  this  problem  is  straight  forward.  The 
data  should  be  high  pass  filtered  before  the  data  is  divided  up 
into  blocks  and  analyzed.  The  frequency  cutoff  should  be  at  1/N 
times  the  sampling  frequency. 

As  a  method  of  implementing  this  filter,  the  following 
method  is  suggested.  First,  the  block  means  are  determined. 
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PERIODIC  REPRESENTATION 


FIGURE  5-4.  LOW  FREQUENCY  BACKGROUND  SUBTRACTION  EXAMPLE 
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Then  values  between  the  blocks  means  at  points  corresponding 
to  the  pixel  positions  are  determined  using  interpolation.  Simple 
linear  interpolation  should  be  sufficient  to  drastically  reduce 
the  problem.  This  interpolated  signal  is  used  as  an  estimate 
of  the  unwanted  low  frequency  components  and  is  subtracted  from 
the  original  signal.  Analysis  of  the  residual  signal  is  then 
performed  as  originally  described. 

Consider  what  happens  when  this  approach  is  applied  to 
the  ramp  signal  example.  Subtraction  of  estimated  low  frequency 
background  results  in  a  null  signal  out.  The  texture  measurements 
will  then  properly  reflect  the  flat  texture  present. 

5. 3  Formal  Definition,  Texture  Matching 

The  texture  matching  procedure  which  has  been  described 
is  formally  summarized  in  this  subsection.  This  is  first  done 
for  the  case  of  a  single  texture  measure.  Next,  this  procedure 
is  generalized  to  handle  multiple  texture  measures.  Again  for 
convenience  one  dimensional  expressions  are  used.- 
5.3.1  Single  Texture  Measure 

A  reference  window  is  selected  from  one  image  to  be 
matched,  and  a  search  area  is  selected  from  the  other  image.  Let 
the  pixel  replacement  resolution  be  given  by  R.  Then  in  pixel 
intensity  space  these  areas  are  represented  by: 

reference  window:  x^,  i  =  1,2,...,R*N 
search  area:  y^,  i  =  1,2,...,R*M,  M>N. 

From  the  pixel  intensity  measurements,  the  texture  measure  is  determined 
for  each  Rth  pixel.  Let  the  texture  measures  derived  for  the  reference 
window  and  the  search  area  be  denoted  by: 

reference  window:  s^,  i  =  1,2,...,N 

t^  ,  '■  ~  1,2,...,M, 


search  area: 


M>N. 


Let  A  be  the  pixel  neighborhood  size  used  to  determine  the  texture 
measure.  Then  each  texture  value  is  determined  using  A  pixel 
values.  Specifically 

si  =  f (x(i-l)*R+l'  x (i-1) *R+2 ' *  *  * '  X(i-1)*R+A) 
ti  “  f  (y(i-l)*R+l'  y (i-1) *R+2 '  *  *  * '  y(i-l)*R+A} 


(5-17) 


(5-18) 


Then  the  normalized  cross  correlation  product  is  formed 
between  the  texture  measures  for  the  reference  window  and  the 
texture  measures  for  the  search  area. 


I  *  lsi  '  W 


j  =  1,2, .. . ,M-N+1 


(5-19) 


with 


3  "  N  i=l  V 
1  N 

h  m  1  i»i  ti+ji 


(5-20) 


(5-21) 


us  -^1  ^  (si'i)2 


(5-22) 


O.  =  1  N  . 


1  «  _  ? 
4  1  (t^-tj* 


'1+3  3 


(5-23) 


The  match  point  is  selected  from  position  of  the 
correlation  peak.  For  good  matches,  the  correlation  peak  value 
will  approach  one.  The  parameters  that  must  be  experimentally 
determined  for  the  matching  of  any  two  types  of  dissimilar  images 
include: 
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a)  the  single  type  of  texture  measure  to  be  used, 

b)  the  size  of  the  reference  window  and  search  area 
to  be  used, 

c)  the  size  of  the  neighborhood  used  to  determine  the 
texture  measure, 

d)  and  the  pixel  replacement  resolution  to  be  used. 

If  a  single  texture  measure  is  to  be  used,  the  local  busyness 
measure  as  defined  by  Equation  4-11  should  be  an  excellent  choice. 
5.3.2  Multiple  Texture  Measures 

To  apply  multiple  texture  measures,  independent  matching 
is  performed  for  each  of  the  measures  used.  Then  an  overall 
correlation  parameter  is  used.  Let  B  give  the  number  of  correlation 
measures  used.  Then  the  overall  correlation  parameter  is  defined  by: 


1  B 
P±  -  |  * 

1  B  u-1 


u 


(5-24) 


The  superscript  u  refers  to  each  of  the  multiple  texture  measures 
used.  The  match  point  is  determined  by  the  peak  of  this  distribution. 
The  normalization  by  B  allows  the  identification  of  good  matches  as 
peak  values  which  approach  one. 

If  it  is  possible  to  determine  that  some  texture  measures 
are  more  important  than  others,  then  a  weighting  procedure  is  easily 
defined.  For  example  the  overall  correlation  parameter  could  be 
defined  by: 
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(5-25) 


„  .  ?  «u  pf, 

pi  “  u=l  1 


with  I  u)U— 1  (5-26) 

u=l 

Here  u>u  is  a  weighting  factor  of  less  than  one  associated  with 
each  of  the  texture  measures.  The  normalization  of  the  weighting 
factors  allows  the  identification  of  good  matches  as  correlation 
peak  values  which  approach  one.  Methods  for  determining  appropriate 
weighting  factors  are  beyond  the  scope  of  this  report. 

5.4  Miscellaneous  Corrections 

In  this  subsection  a  variety  of  special  problems  which 
can  interfere  with  texture  matching  are  discussed.  Methods  are 
discussed  for  dealing  with  data  noise  characteristics,  different 
detector  spatial  acceptances,  drastically  different  intensity 
distributions,  and  coherency  effects  found  with  imaging  radars. 

Detailed  algorithms  for  dealing  with  these  problems  are  not,  in 
general,  presented.  The  emphasis  here  is  to  point  out  the  potential 
existance  of  these  problems  and  to  suggest  lines  along  which  the 
problems  could  be  dealt  with. 

5.4.1  Noise  Characteristics 

Image  data  is  frequently  corrupted  by  significant  amount 
of  noise.  This  noise  can  interfere  with  the  image  matching  process. 

We  will  consider  here  only  the  most  commonly  found  type  of  noise, 
that  of  white,  Gaussian  noise.  Consider  an  example  in  which 
the  two  images  to  be  matched  have  differing  levels  of  Gaussian 
noise.  Assume  that  each  intensity  measurement  from  one  image  has 
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a  noise  standard  deviation  component  which  is  twice  as  large  as  that 
in  the  second  image.  Consider  what  happens  in  a  region  of  the  scene 
which  is  of  constant  intensity  if  not  corrupted  by  noise.  Using  the 
texture  matching  approach,  the  texture  measures  derived  from  spectral 
analysis  from  this  region  for  the  first  image  will  be  on  the  average 
twice  as  large  as  texture  measures  from  this  region  in  the  second 
image.  For  a  match  attempt  with  a  reference  window  which  includes 
this  region,  the  match  will  be  adversely  affected  because  of  this. 

A  straight  forward  solution  to  this  problem  exists.  Noise  should 
be  added  to  the  second  image  to  bring  it  up  to  the  level  of  the  first 
image.  This  paradoxial  solution  gives  texture  measures  in  the  region 
in  question  which  are  now  on  the  average  of  equal  magnitude. 

5.4.2  Spatial  Acceptance 

Differences  in  spatial  acceptances  of  the  dissimilar 
images  to  be  matched  will  interfere  with  the  determination  of 
corresponding  points  using  texture  matching.  Differences  in 
spatial  acceptance  causes  each  of  the  images  to  be  blurred  by 
different  amounts.  This  causes  the  texture  at  corresponding  points 
to  be  different  from  one  image  to  the  next.  Again  the  solution  to 
this  problem  is  fairly  straight  forward. 

For  simplicity  consider  the  following  example.  A  single 
aerial  camera  has  been  used  to  produce  two  images  of  the  same 
area.  The  only  difference  between  the  two  images  is  the  altitude 
used.  The  second  image  io  taken  from  an  altitude  which  is  twice 
as  high  as  the  altitude  used  for  the  first  image.  After  sampling 
both  images  with  what  correspondes  to  equal  sample  spacing  on 
the  ground,  the  two  images  will  not  be  quite  identical.  The  second 
image  will  be  blurred  with  respect  to  the  first.  In  practice  it 
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should  still  be  possible  to  obtain  match  points  with  methods  that 
have  been  described.  However,  it  is  possible  to  improve  the  matching 
process . 

Improved  results  are  obtained  by  adding  additional 
blurring  to  the  first  image.  The  blurring  can  be  accomplished 
with  the  use  of  an  appropriate  low  pass  filter.  The  required 
properties  of  the  low  pass  filter,  can  be  easily  determined  if 
the  spatial  acceptances  for  the  two  images  are  known.  This  is 
illustrated  for  the  example  being  considered  in  Figure  5-5.  The 
spatial  acceptances  for  the  two  images  are  shown  in  the  form  of 
their  modulation  transfer  functions.  A  modulation  transfer  function 
gives  the  attenuation  of  spatial  waves  as  a  function  of  frequency. 
Frequencies  have  been  expressed  by  normalizing  by  the  sampling 
frequency.  The  frequency  response  of  the  required  filter  is  given 
by  the  ratio  of  the  acceptance  of  second  image  to  the  acceptance 
of  the  first  image  as  a  function  of  frequency.  The  required 
response  for  the  example  is  shown  in  Figure  5-5. 

5.4.3  Intensity  Distributions 

In  general,  the  matching  procedure  that  has  been  described 
should  be  quite  insensitive  to  broad  differences  in  the  magnitude  and  the 
shape  of  the  intensity  distributions  of  the  two  images  to  be  matched 
However,  with  the  matching  of  an  optical  image  with  a  radar  image 
extreme  differences  can  be  found  in  their  intensity  distributions. 

In  general,  the  intensity  distribution  of  an  optical  image  is 
relatively  flat.  On  the  other  hand,  the  intensity  distribution  of 
a  typical  radar  image  consists  of  an  extremely  large  clutter  peak 
at  /ery  low  intensities  plus  tail  made  up  of  the  intensities  of  a 
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FIGURE  5-5.  EXAMPLE  MATCHING  SPATIAL  ACCEPTANCE 
WITH  LOW  PASS  FILTERING. 


few  very  bright  targets.  Well  over  half  of  the  intensity  values 
will  frequently  fall  in  the  clutter  peak.  Typical  distributions 
are  illustrated  in  Figure  5-6. 

Under  these  circumstances  it  is  suggested  that  improved 
matching  results  could  be  obtained  by  modifying  the  radar  image 
intensity  distribution  with  the  use  of  an  intensity  transformation. 

The  transform  could  be  of  the  form 

1'  =  alb  (5-27) 

where  b  is  less  than  one.  Figure  5-6  shows  a  modified  intensity 
distribution  of  the  original  radar  distribution  that  could  be 
obtained  by  a  transformation  of  this  form. 

5.4.4  Coherency  Effects 

The  matching  of  high  resolution  radar  images  gives  a 
severe  test  to  any  image  matching  technique.  The  coherent  nature 
of  imaging  radar  produces  targets  with  a  high  frequency  scintillation 
component.  These  waves  shift  in  position  with  the  slightest  change 
in  taking  parameters .  Because  of  this ,  conventional  correlation 
matching  procedures  will  frequently  not  work  on  radar  images  taken 
with  nominally  the  same  taking  parameters.  For  this  particular 
case,  the  texture  matching  procedure  should  overcome  this  difficulty. 
However,  between  imaging  radars  of  different  resolutions  or  between 
a  radar  image  and  an  optical  image,  coherent  effects  present  a 
difficult  problem. 

One  suggestion  to  deal  with  this  problem  is  to  treat  the 
high  frequency  coherent  components  as  noise.  Given  that  the  high 
frequency  components  of  one  of  the  images  are  dominated  by  noise, 
matching  using  texture  matching  will  be  improved  by  low  pass  filtering 
both  the  images.  The  low  pass  filter  should  cutoff  at  the  frequency 

at  which  noise  starts  to  dominate. 

92 


5.5 


Match  Point  Accuracy 

In  general,  match  point  accuracy  will  vary  as  match  points 
are  obtained  across  a  pair  of  images.  In  some  areas  of  a  scene  it 
may  not  be  possible  to  obtain  well  defined  match  points,  while  in 
other  areas  highly  accurate  match  points  will  be  obtained.  There 
are  several  approaches  which  can  be  taken  to  predict  or  evaluate 
the  accuracy  of  a  possible  match  point.  In  general,  the  more  active 
a  local  area  of  the  scene  is,  the  more  accurate  the  match  point 
determinations  will  be.  In  an  area  completely  devoid  of  any 
activity,  match  point  determination  will  not  be  possible.  In  an 
area  dominated  by  sharp  edges  or  other  sharply  defined  structures, 
highly  accurate  match  points  will  be  obtained.  This  fact  can  be 
used  to  select  areas  of  the  scene  which  give  the  most  accurate  match 
points.  Local  activity  could  be  measured  by  using  the  local 
busyness  texture  measure  (tg) .  This  approach  could  be  very  effectively 
used  in  a  system  designed  to  approximately  register  images  with  the 
measurement  of  a  few  match  points. 

Once  a  match  has  been  undertaken,  information  generated 
in  the  matching  process  can  be  used  to  further  refine  the  estimate 
of  match  point  accuracy.  The  value  of  the  correlation  parameter  at 
the  correlation  peak  gives  a  measure  of  texture  similarity  at  the 
match  point.  A  threshold  can  be  set  to  require  a  certain  level  of 
texture  similarity  for  a  good  match.  The  accuracy  of  match  can  be 
evaluated  using  the  sharpness  of  the  correlation  peak.  In  general, 
the  sharper  the  correlation  peak,  the  more  accurate  the  match. 


94 


6.0 


SUMMARY 


To  solve  the  general  matching  problem  of  dissimilar  images, 
a  texture  matching  approach  has  been  taken.  It  is  assumed  that  from 
one  image  of  a  scene  to  another  there  are  elements  of  texture  which 
will  remain  invariant.  This  approach  can  potentially  be  used  in  a 
wide  range  of  image  matching  applications.  These  range  from  the 
approximate  registration  problem  where  a  few  match  points  are 
required  to  the  problem  of  determining  thousands  of  match  points. 

The  key  to  this  approach  is  the  development  of  specific  conversion 
algorithms  for  the  pixels  of  each  system  of  interest,  so  that  the 
resultant  pixel  values  are  effectively  system  independent.  The 
actual  pixel  replacement  algorithms  will  have  to  be  empirically 
tailored  for  each  specific  collection  system. 

To  complement  this  completely  automatic  approach,  a 
semi-automatic  approach  is  presented  which  can  be  used  for  the 
approximate  registration  problem.  This  approach,  computer  assisted 
manual  image  registration  (CAMIR) ,  provides  computer  and  special- 
purpose  hardware  assistance  to  the  operator  who  is  required  to 
measure  matching  points  from  each  of  the  images.  This  system 
could  be  readily  constructed  given  the  current  state-of-the-art 
of  digital  hardware. 

The  first  step  in  the  match  problem  is  to  apply  geometric 
corrections  using  the  initial  knowledge  of  the  taking  parameters. 

The  inclusion  of  geometric  effects  requires  three  steps: 
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a)  A  grid  of  points  is  defined  on  the  ground  at  positions 
at  which  it  is  desired  to  know  the  intensity  values, 

b)  These  grid  points  are  mapped  into  the  projection 
space  and  their  corresponding  positions  in  projection 
space  are  determined. 

c)  Pixel  intensity  values  at  these  positions  in  pro¬ 
jection  space  are  determined  using  interpolation  and 
the  existing  intensity  measurements. 

The  geometrically  corrected  image  is  formed  from  these  interpolated 
pixel  values.  After  geometric  correction,  the  pair  of  images  to 
be  matched  will  usually  differ  by  a  translational  shift.  This 
misregistration  is  due  to  residual  uncertainties  in  the  aim 
point  determination.  Rotational  and  scale  differences  between 
the  images  to  be  matched  will  normally  be  small.  Application 
of  geometric  corrections  greatly  reduces  the  dimensions  of  the 
matching  problem. 

The  CAMIR  system  makes  use  of  the  exceptional  ability 
of  a  human  operator  to  quickly  solve  complex  pattern  recognition 
problems.  Computer  and  special  high-speed  digital  hardware  are 
combined  in  a  manner  to  permit  an  operator  to  solve  the  registra¬ 
tion  problem  with  maximum  efficiency.  The  operator  has  a  CRT 
display  to  view  each  of  the  two  dissimilar  images.  His  problem 

is  to  find  and  measure  three  or  more  match  points  from  the  extremes 
of  the  images.  The  hardware  permits  the  operator  to  have  complete 
freedom  to  explore  and  examine  the  images.  He  can  translate, 
rotate,  and  magnify  or  demagnify  the  sub-image  being  viewed. 

The  ability  to  control  image  sharpness  and  noise  can  be  provided 
under  the  real-time  control  of  the  operator.  This  is  provided 
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by  adjusting  spatial  response  with  a  modulation  transfer  function 
compensator  (MTFC) .  Similarly,  the  ability  to  modify  scene 
intensities  can  be  provided  under  the  real-time  control  of  the 
operator.  This  dynamic  range  adjustment  (DRA)  maps  intensity 
levels  in  any  desired  manner. 

To  provide  for  maximum  efficiency  in  registering  images, 
geometric  corrections  are  always  applied  before  the  images  are 
displayed.  To  be  able  to  apply  corrections  to  a  variety  of  taking 
systems,  the  transformations  from  object  space  to  projection  space 
are  expressed  in  a  general  polynomial  form.  With  this  approach 
software  can  cast  the  correction  problem  for  any  system  into  a 
uniform  form.  Special  purpose  hardware  can  then  solve  the  problem 
in  real  time.  Once  the  three  corresponding  points  have  been 
measured  by  the  operator,  the  images  are  brought  together  into 
registration.  This  is  accomplished  by  adjusting  the  parameters 
of  the  distortion  polynomials. 

This  system  could  be  used  as  an  integral  part  of  a 
variety  of  systems  designed  to  attack  map  making  problems  with 
digital  techniques.  These  could  include  map  updating  systems 
or  systems  designed  to  build  up  contour  data  base  information. 

This  system  can  be  built  using  currently  available  digital 
components.  For  E-Systems  to  build  such  a  system,  no  research  or 
development  effort  would  be  required.  The  only  requirement  to  build 
the  image  manipulator  would  be  engineering  effort. 
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The  texture  matching  approach  to  image  matching  has  the 
advantage  of  being  completely  automatic.  The  first  step  in  the 
process  is  to  correct  for  all  known  geometric  distortions. 

After  geometric  correction,  the  two  images  to  be  matched  will 
normally  differ  in  registration  by  a  translation  shift.  In  the 
next  step,  each  of  the  images  to  be  matched  is  transformed  into 
a  texture  measure  space.  Each  pixel  value  is  replaced  by  some 
local  measure  of  the  surround  (texture) .  It  is  not  always 
necessary  to  replace  each  and  every  pixel  intensity  value  by  a 
texture  measure,  and  the  frequency  of  replacement  that  is  used 
is  referred  to  as  the  pixel  replacement  resolution.  The  match 
point  is  determined  by  correlating  in  the  texture  measure  space. 

A  reference  window  is  chosen  from  one  image  and  a  larger  search 
area  is  chosen  from  the  other  image.  The  normalized  cross 
correlation  product  is  used  and  the  match  point  is  determined  by 
the  position  of  the  correlation  peak.  A  threshold  can  be  set 
using  the  value  of  the  correlation  parameter  at  the  peak  to 
determine  if  a  good  match  has  been  found.  This  process  can 
be  repeated  across  the  two  images  to  obtain  match  points  as 
frequently  as  desired.  This  technique  is  generalized  to  permit 
the  use  of  multiple  texture  measures . 

Local  spectral  analysis  is  an  extremely  effective  tool 
that  can  be  used  to  develop  appropriate  texture  measures.  The 
local  spectral  analysis  can  be  performed  with  the  use  of  the  two 
dimensional  discrete  Fourier  transform.  This  transform  has  the 
advantage  of  having  a  fast  algorithm  (the  FFT)  available  for 


its  implementation.  Texture  measures  can  be  defined  by  grouping 
together  similar  spectral  components  using  the  general  prescription: 
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Here  Xu  v  represents  th  complex  spectral  components.  Spectral 
components  can  be  grouped  according  to  their  radial  frequency 
(  n|u2+v2) .  A  variety  of  possible  texture  measures  are  given  for 
the  four-by-four  and  eight-by-eight  neighborhood  sizes.  The  size 
of  the  neighborhood  defines  the  local  area  used  in  determining 
the  texture  measures. 

If  it  is  desired  to  use  a  single  measure  of  texture, 
the  local  busyness  measure  should  be  extremely  effective.  This 
is  the  texture  measure  which  is  used  to  demonstrate  the  texture 
matching  approach.  It  is  defined  as  the  square  root  of  the  AC 
energy  present  in  the  pixel  neighborhood. 

To  use  texture  matching,  the  algorithm  must  be  empirically 
tailored  to  the  combination  of  collection  systems  that  are  used. 

The  parameters  which  must  be  experimentally  determined  include: 

a)  the  selection  of  the  texture  measure  or  measures 
to  be  used, 

b)  the  size  of  the  reference  window  and  search  area  to 
be  used, 

c)  the  size  of  the  neighborhood  used  to  determine  the 
texture  measure, 

d)  and  the  pixel  replacement  resolution  to  be  used. 
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There  are  a  variety  of  special  problems  or  effects 
which  can  interfere  with  the  matching  process.  A  number  of  these 
effects  are  discussed  and  methods  of  dealing  with  these  problems 
are  presented.  These  problems  include  data  noise  characteristics, 
different  detector  spatial  acceptances,  drastically  different 
intensity  distributions,  and  coherency  effects  found  in  dealing 
with  imaging  radars. 


